Create a new analysis directory...
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] "/Users/swvanderlaan/git/CirculatoryHealth/AE_20211201_YAW_SWVANDERLAAN_HDAC9"
[1] "_archived" "1. AEDB.CEA.baseline.nb.html" "1. AEDB.CEA.baseline.Rmd"
[4] "2. SNP_analyses.Rmd" "20220319.HDAC9.AEDB.CEA.baseline.RData" "20220319.HDAC9.AESCRNA.results.RData"
[7] "20220319.HDAC9.bulkRNAseq.main_analysis.RData" "20220319.HDAC9.bulkRNAseq.preparation.RData" "3.1 bulkRNAseq.preparation.nb.html"
[10] "3.1 bulkRNAseq.preparation.Rmd" "3.2 bulkRNAseq.main_analysis.nb.html" "3.2 bulkRNAseq.main_analysis.Rmd"
[13] "3.3 bulkRNAseq.additional_figures.Rmd" "4. scRNAseq.nb.html" "4. scRNAseq.Rmd"
[16] "AE_20211201_YAW_SWVANDERLAAN_HDAC9.Rproj" "AnalysisPlan" "HDAC9"
[19] "images" "LICENSE" "README.html"
[22] "README.md" "references.bib" "renv"
[25] "renv.lock" "scripts" "SNP"
[28] "targets" source(paste0(PROJECT_loc, "/scripts/functions.R"))install.packages.auto("pander")
install.packages.auto("readr")
install.packages.auto("optparse")
install.packages.auto("tools")
install.packages.auto("dplyr")
install.packages.auto("tidyr")
install.packages.auto("naniar")
# To get 'data.table' with 'fwrite' to be able to directly write gzipped-files
# Ref: https://stackoverflow.com/questions/42788401/is-possible-to-use-fwrite-from-data-table-with-gzfile
# install.packages("data.table", repos = "https://Rdatatable.gitlab.io/data.table")
library(data.table)
install.packages.auto("tidyverse")
install.packages.auto("knitr")
install.packages.auto("DT")
install.packages.auto("eeptools")
install.packages.auto("openxlsx")
install.packages.auto("haven")
install.packages.auto("tableone")
install.packages.auto("sjPlot")
install.packages.auto("BlandAltmanLeh")
# Install the devtools package from Hadley Wickham
install.packages.auto('devtools')
# for plotting
install.packages.auto("pheatmap")
install.packages.auto("forestplot")
install.packages.auto("ggplot2")
install.packages.auto("ggpubr")
install.packages.auto("UpSetR")
devtools::install_github("thomasp85/patchwork")Using github PAT from envvar GITHUB_PAT
Skipping install of 'patchwork' from a github remote, the SHA1 (79223d30) has not changed since last install.
Use `force = TRUE` to force installationinstall.packages.auto("GGally")
install.packages.auto("ggcorrplot")
# for Seurat etc
install.packages.auto("GenomicFeatures")
install.packages.auto("GenomicRanges")
install.packages.auto("SummarizedExperiment")
install.packages.auto("DESeq2")
install.packages.auto("org.Hs.eg.db")
install.packages.auto("mygene")
install.packages.auto("TxDb.Hsapiens.UCSC.hg19.knownGene")
install.packages.auto("org.Hs.eg.db")
install.packages.auto("AnnotationDbi")
install.packages.auto("EnsDb.Hsapiens.v86")Loading required package: EnsDb.Hsapiens.v86
Loading required package: ensembldb
Attaching package: 'ensembldb'
The following object is masked from 'package:openxlsx':
addFilter
The following object is masked from 'package:dplyr':
filter
The following object is masked from 'package:stats':
filterinstall.packages.auto("EnhancedVolcano")
Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")
### UtrechtScienceParkColoursScheme
###
### WebsitetoconvertHEXtoRGB:http://hex.colorrrs.com.
### Forsomefunctionsyoushoulddividethesenumbersby255.
###
### No. Color HEX (RGB) CHR MAF/INFO
###---------------------------------------------------------------------------------------
### 1 yellow #FBB820 (251,184,32) => 1 or 1.0>INFO
### 2 gold #F59D10 (245,157,16) => 2
### 3 salmon #E55738 (229,87,56) => 3 or 0.05<MAF<0.2 or 0.4<INFO<0.6
### 4 darkpink #DB003F ((219,0,63) => 4
### 5 lightpink #E35493 (227,84,147) => 5 or 0.8<INFO<1.0
### 6 pink #D5267B (213,38,123) => 6
### 7 hardpink #CC0071 (204,0,113) => 7
### 8 lightpurple #A8448A (168,68,138) => 8
### 9 purple #9A3480 (154,52,128) => 9
### 10 lavendel #8D5B9A (141,91,154) => 10
### 11 bluepurple #705296 (112,82,150) => 11
### 12 purpleblue #686AA9 (104,106,169) => 12
### 13 lightpurpleblue #6173AD (97,115,173/101,120,180) => 13
### 14 seablue #4C81BF (76,129,191) => 14
### 15 skyblue #2F8BC9 (47,139,201) => 15
### 16 azurblue #1290D9 (18,144,217) => 16 or 0.01<MAF<0.05 or 0.2<INFO<0.4
### 17 lightazurblue #1396D8 (19,150,216) => 17
### 18 greenblue #15A6C1 (21,166,193) => 18
### 19 seaweedgreen #5EB17F (94,177,127) => 19
### 20 yellowgreen #86B833 (134,184,51) => 20
### 21 lightmossgreen #C5D220 (197,210,32) => 21
### 22 mossgreen #9FC228 (159,194,40) => 22 or MAF>0.20 or 0.6<INFO<0.8
### 23 lightgreen #78B113 (120,177,19) => 23/X
### 24 green #49A01D (73,160,29) => 24/Y
### 25 grey #595A5C (89,90,92) => 25/XY or MAF<0.01 or 0.0<INFO<0.2
### 26 lightgrey #A2A3A4 (162,163,164) => 26/MT
###
### ADDITIONAL COLORS
### 27 midgrey #D7D8D7
### 28 verylightgrey #ECECEC"
### 29 white #FFFFFF
### 30 black #000000
###----------------------------------------------------------------------------------------------
uithof_color = c("#FBB820","#F59D10","#E55738","#DB003F","#E35493","#D5267B",
"#CC0071","#A8448A","#9A3480","#8D5B9A","#705296","#686AA9",
"#6173AD","#4C81BF","#2F8BC9","#1290D9","#1396D8","#15A6C1",
"#5EB17F","#86B833","#C5D220","#9FC228","#78B113","#49A01D",
"#595A5C","#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")
uithof_color_legend = c("#FBB820", "#F59D10", "#E55738", "#DB003F", "#E35493",
"#D5267B", "#CC0071", "#A8448A", "#9A3480", "#8D5B9A",
"#705296", "#686AA9", "#6173AD", "#4C81BF", "#2F8BC9",
"#1290D9", "#1396D8", "#15A6C1", "#5EB17F", "#86B833",
"#C5D220", "#9FC228", "#78B113", "#49A01D", "#595A5C",
"#A2A3A4", "#D7D8D7", "#ECECEC", "#FFFFFF", "#000000")
### ----------------------------------------------------------------------------Here we perform the main analyses. We first load the prepared RNAseq data, and extract only the relevant genes.
Here we obtain data from the HDAC9 in plaques.
library(openxlsx)
gene_list_df <- read.xlsx(paste0(PROJECT_loc, "/targets/Genes.xlsx"), sheet = "Genes")
target_genes <- unlist(gene_list_df$Gene)
target_genes[1] "HDAC9" "TWIST1" "IL6" "IL1B" We simply load the previously saved RDS-file and extract
the clinical and RNAseq data from that.
AEDB.CEA <- readRDS(file = paste0(OUT_loc, "/20220319.",TRAIT_OF_INTEREST,".AEDB.CEA.RDS"))
AEDB.CEA$STUDY_NUMBER <- paste0("ae", AEDB.CEA$STUDY_NUMBER)
head(AEDB.CEA$STUDY_NUMBER)[1] "ae3" "ae7" "ae9" "ae12" "ae13" "ae15"AERNASE <- readRDS(file = paste0(OUT_loc, "/20220319.AERNA.CEA.623pts.SE.after_qc.IC_academic.RDS"))AERNASE.clin <- as_tibble(colData(AERNASE))
AERNASE.hdac9 <- subset(AERNASE, subset = (rowRanges(AERNASE)$symbol %in% target_genes))
row.names(AERNASE.hdac9) <- rowData(AERNASE.hdac9)$symbol
temp <- as_tibble(t(assay(AERNASE.hdac9)), rownames = "STUDY_NUMBER")
AERNASE.clin.hdac9 <- as.data.frame(merge(AERNASE.clin, temp, by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE))
rm(temp)Getting some summary statistics.
cat("Average expression of a random selection of 1000 genes.\n")Average expression of a random selection of 1000 genes.set.seed(141619)
mean(rowMeans(as_tibble(assay(AERNASE)) %>% dplyr::sample_n(1000)))[1] 12.54526cat("\nAverage expression of target genes.\n")
Average expression of target genes.rowMeans(assay(AERNASE.hdac9)) HDAC9 TWIST1 IL1B IL6
28.3226324 5.5553772 5.6950241 0.4109149 cat("\nGene information.\n")
Gene information.rowRanges(AERNASE.hdac9)GRanges object with 4 ranges and 3 metadata columns:
seqnames ranges strand | feature_id symbol GENEBIOTYPE_EnsDb86
<Rle> <IRanges> <Rle> | <character> <character> <character>
HDAC9 7 18086949-19002416 + | ENSG00000048052 HDAC9 protein_coding
TWIST1 7 19020991-19117636 - | ENSG00000122691 TWIST1 protein_coding
IL1B 2 112829751-112836816 - | ENSG00000125538 IL1B protein_coding
IL6 7 22725884-22732002 + | ENSG00000136244 IL6 protein_coding
-------
seqinfo: 338 sequences from an unspecified genome; no seqlengthsAERNASE.hdac9.genedata <- as_tibble(rowRanges(AERNASE.hdac9))
AERNASE.hdac9.genedataThe analyses are focused on three elements:
Age]Gender]Hypertension.composite]DiabetesStatus]SmokerStatus]LDL_final]Med.Statin.LLD]Med.all.antiplatelet]GFR_MDRD]BMI]MedHx_CVD] combination of
[CAD_history, Stroke_history,
Peripheral.interv]stenose]ORdate_year] as we discovered in Van
Lammeren et al. the composition of the plaque and therefore the
Athero-Express Biobank Study has changed over the years. Likely through
changes in lifestyle and primary prevention regimes.We will analyze the data through four different models
In the cross-sectional analysis of plaque MCP1 levels we will focus on the following plaque vulnerability phenotypes:
We inspect the plaque characteristics, and
inverse-rank normal transformation continuous
phenotypes.
# macrophages
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$MAC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.0937 -0.3666 0.3275 0.2830 0.9611 2.6777 15 min_macmean <- min(AERNASE.clin.hdac9$MAC_rankNorm, na.rm = TRUE)
cat(paste0("\nMinimum value % macrophages: ",min_macmean,".\n"))
Minimum value % macrophages: -2.09371245439672.ggpubr::gghistogram(AERNASE.clin.hdac9, "MAC_rankNorm",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "% macrophages",
xlab = "inverse-rank normalized %",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 15 rows containing non-finite values (stat_bin).
# smooth muscle cells
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$SMC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.64372 -0.65587 -0.02105 -0.03204 0.63927 2.55750 17 min_smcmean <- min(AERNASE.clin.hdac9$SMC_rankNorm, na.rm = TRUE)
cat(paste0("\nMinimum value % smooth muscle cells: ",min_smcmean,".\n"))
Minimum value % smooth muscle cells: -2.64372511353808.ggpubr::gghistogram(AERNASE.clin.hdac9, "SMC_rankNorm",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "% smooth muscle cells",
xlab = "inverse-rank normalized %",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 17 rows containing non-finite values (stat_bin).
# vessel density
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$VesselDensity_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-1.9203 -0.4576 0.1109 0.1365 0.7694 3.1703 53 min_vesseldensity <- min(AERNASE.clin.hdac9$VesselDensity_rankNorm, na.rm = TRUE)
min_vesseldensity[1] -1.920282cat(paste0("\nMinimum value number of intraplaque neovessels per 3-4 hotspots: ",min_vesseldensity,".\n"))
Minimum value number of intraplaque neovessels per 3-4 hotspots: -1.92028163343836.ggpubr::gghistogram(AERNASE.clin.hdac9, "VesselDensity_rankNorm",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "number of intraplaque neovessels per 3-4 hotspots",
xlab = "inverse-rank normalized number",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 53 rows containing non-finite values (stat_bin).
Given their strong correlation, we also introduce a macrophages/smooth muscle cell ratio. This is a proxy of the extend to which a plaque is inflammed (‘unstable’) as compared to ‘stable’.
AERNASE.clin.hdac9$MAC_SMC_ratio <- AERNASE.clin.hdac9$MAC_rankNorm / AERNASE.clin.hdac9$SMC_rankNorm
AERNASE.clin.hdac9$MAC_SMC_ratio_rank <- qnorm((rank(AERNASE.clin.hdac9$MAC_SMC_ratio, na.last = "keep") - 0.5) / sum(!is.na(AERNASE.clin.hdac9$MAC_SMC_ratio)))
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$MAC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.0937 -0.3666 0.3275 0.2830 0.9611 2.6777 15 summary(AERNASE.clin.hdac9$SMC_rankNorm) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-2.64372 -0.65587 -0.02105 -0.03204 0.63927 2.55750 17 summary(AERNASE.clin.hdac9$MAC_SMC_ratio_rank) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-3.1469 -0.6732 0.0000 0.0000 0.6732 3.1469 17 ggpubr::gghistogram(AERNASE.clin.hdac9, "MAC_SMC_ratio_rank",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "median",
#add_density = TRUE,
rug = TRUE,
#add.params = list(color = "black", linetype = 2),
title = "macrophages/smooth muscle cells ratio",
xlab = "inverse-rank normalized",
ggtheme = theme_minimal())Warning: Using `bins = 30` by default. Pick better value with the argument `bins`.
Warning: Removed 17 rows containing non-finite values (stat_bin).
library(dplyr)
# calcification
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Calc.bin) no/minor moderate/heavy NA's
285 330 8 contrasts(AERNASE.clin.hdac9$Calc.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$CalcificationPlaque <- as.factor(AERNASE.clin.hdac9$Calc.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(CalcificationPlaque)) %>%
dplyr::group_by(Gender, CalcificationPlaque) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "CalcificationPlaque", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Calcification",
xlab = "calcification",
ggtheme = theme_minimal())
rm(df)
# collagen
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Collagen.bin) no/minor moderate/heavy NA's
118 496 9 contrasts(AERNASE.clin.hdac9$Collagen.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$CollagenPlaque <- as.factor(AERNASE.clin.hdac9$Collagen.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(CollagenPlaque)) %>%
dplyr::group_by(Gender, CollagenPlaque) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "CollagenPlaque", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Collagen",
xlab = "collagen",
ggtheme = theme_minimal())
rm(df)
# fat 10%
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Fat.bin_10) <10% >10% NA's
145 470 8 contrasts(AERNASE.clin.hdac9$Fat.bin_10) >10%
<10% 0
>10% 1AERNASE.clin.hdac9$Fat10Perc <- as.factor(AERNASE.clin.hdac9$Fat.bin_10)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(Fat10Perc)) %>%
dplyr::group_by(Gender, Fat10Perc) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "Fat10Perc", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Intraplaque fat",
xlab = "intraplaque fat",
ggtheme = theme_minimal())
rm(df)
# macrophages binned
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$Macrophages.bin) no/minor moderate/heavy NA's
265 347 11 contrasts(AERNASE.clin.hdac9$Macrophages.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$MAC_binned <- as.factor(AERNASE.clin.hdac9$Macrophages.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(MAC_binned)) %>%
dplyr::group_by(Gender, MAC_binned) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "MAC_binned", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Macrophages (binned)",
xlab = "Macrophages",
ggtheme = theme_minimal())
rm(df)
# SMC binned
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$SMC.bin) no/minor moderate/heavy NA's
194 418 11 contrasts(AERNASE.clin.hdac9$SMC.bin) moderate/heavy
no/minor 0
moderate/heavy 1AERNASE.clin.hdac9$SMC_binned <- as.factor(AERNASE.clin.hdac9$SMC.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(SMC_binned)) %>%
dplyr::group_by(Gender, SMC_binned) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "SMC_binned", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "SMC (binned)",
xlab = "SMC",
ggtheme = theme_minimal())
rm(df)
# IPH
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$IPH.bin) no yes NA's
235 379 9 contrasts(AERNASE.clin.hdac9$IPH.bin) yes
no 0
yes 1AERNASE.clin.hdac9$IPH <- as.factor(AERNASE.clin.hdac9$IPH.bin)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(IPH)) %>%
dplyr::group_by(Gender, IPH) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "IPH", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Intraplaque hemorrhage",
xlab = "intraplaque hemorrhage",
ggtheme = theme_minimal())
rm(df)
# Symptoms
cat("Summary of data.\n")Summary of data.summary(AERNASE.clin.hdac9$AsymptSympt) Asymptomatic Ocular and others Symptomatic
80 139 404 contrasts(AERNASE.clin.hdac9$AsymptSympt) Ocular and others Symptomatic
Asymptomatic 0 0
Ocular and others 1 0
Symptomatic 0 1AERNASE.clin.hdac9$AsymptSympt <- as.factor(AERNASE.clin.hdac9$AsymptSympt)
df <- AERNASE.clin.hdac9 %>%
dplyr::filter(!is.na(AsymptSympt)) %>%
dplyr::group_by(Gender, AsymptSympt) %>%
dplyr::summarise(counts = n()) `summarise()` has grouped output by 'Gender'. You can override using the `.groups` argument.ggpubr::ggbarplot(df, x = "AsymptSympt", y = "counts",
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#DB003F", "#1290D9"),
label = TRUE, lab.vjust = 2, lab.col = "#FFFFFF",
title = "Symptoms",
xlab = "symptoms",
ggtheme = theme_minimal())
rm(df)Here we compare the HDAC9 plaque gene expression levels.
p1 <- ggpubr::ggscatter(AERNASE.clin.hdac9,
x = "ORyear",
y = "HDAC9",
color = "#1290D9",
# fill = "Gender",
# palette = c("#1290D9", "#DB003F"),
add = "reg.line",
add.params = list(color = "black", linetype = 2),
cor.coef = TRUE, cor.method = "spearman",
xlab = "year of surgery",
ylab = "experiment 1",
title = "HDAC9 (normalized expression)",
ggtheme = theme_minimal())
p1 `geom_smooth()` using formula 'y ~ x'
rm(p1)In this section we make some variables to assist with analysis.
We fix the diabetes status variable.
# Fix diabetes
attach(AEDB.CEA)
AEDB.CEA[,"DiabetesStatus"] <- NA
AEDB.CEA$DiabetesStatus[DM.composite == -999] <- NA
AEDB.CEA$DiabetesStatus[DM.composite == 0] <- "Control (no Diabetes Dx/Med)"
AEDB.CEA$DiabetesStatus[DM.composite == 1] <- "Diabetes"
detach(AEDB.CEA)
table(AEDB.CEA$DM.composite, AEDB.CEA$DiabetesStatus)
Control (no Diabetes Dx/Med) Diabetes
0 1820 0
1 0 574# AEDB.CEA.temp <- subset(AEDB.CEA, select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "DM.composite", "DiabetesStatus"))
# require(labelled)
# AEDB.CEA.temp$Gender <- to_factor(AEDB.CEA.temp$Gender)
# AEDB.CEA.temp$Hospital <- to_factor(AEDB.CEA.temp$Hospital)
# AEDB.CEA.temp$Artery_summary <- to_factor(AEDB.CEA.temp$Artery_summary)
# AEDB.CEA.temp$DiabetesStatus <- to_factor(AEDB.CEA.temp$DiabetesStatus)
#
# DT::datatable(AEDB.CEA.temp[1:10,], caption = "Excerpt of the whole AEDB.CEA.", rownames = FALSE)
#
# rm(AEDB.CEA.temp)We will also fix a history of CAD, stroke or peripheral intervention
status variable. This will be based on CAD_history,
Stroke_history, and Peripheral.interv
# AEDB.CEA$CAD_history
# AEDB.CEA$Stroke_history
# AEDB.CEA$Peripheral.interv
# Fix diabetes
attach(AEDB.CEA)
AEDB.CEA[,"MedHx_CVD"] <- NA
AEDB.CEA$MedHx_CVD[CAD_history == "No history CAD" | Stroke_history == 0 | Peripheral.interv == "no"] <- "No"
AEDB.CEA$MedHx_CVD[CAD_history == "History CAD" | Stroke_history == 1 | Peripheral.interv == "yes"] <- "yes"
detach(AEDB.CEA)
table(AEDB.CEA$CAD_history)
Missing No history CAD History CAD
0 1618 756 table(AEDB.CEA$Stroke_history)
0 1
1642 779 table(AEDB.CEA$Peripheral.interv)
no yes
1868 482 table(AEDB.CEA$MedHx_CVD)
No yes
892 1529 # AEDB.CEA.temp <- subset(AEDB.CEA, select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "diet810", "AlcoholUse"))
# require(labelled)
# AEDB.CEA.temp$Gender <- to_factor(AEDB.CEA.temp$Gender)
# AEDB.CEA.temp$Hospital <- to_factor(AEDB.CEA.temp$Hospital)
# AEDB.CEA.temp$Artery_summary <- to_factor(AEDB.CEA.temp$Artery_summary)
# AEDB.CEA.temp$AlcoholUse <- to_factor(AEDB.CEA.temp$AlcoholUse)
#
# DT::datatable(AEDB.CEA.temp[1:10,], caption = "Excerpt of the whole AEDB.CEA.", rownames = FALSE)
#
# rm(AEDB.CEA.temp)AERNASE.clin.hdac9$DiabetesStatus <- NULL
AERNASE.clin.hdac9$MedHx_CVD <- NULL
AERNASE.clin.hdac9 <- merge(AERNASE.clin.hdac9,
subset(AEDB.CEA, select = c("STUDY_NUMBER", "dateok", "DiabetesStatus", "MedHx_CVD")),
by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE)
AERNASE.clin.hdac9$MedHx_CVD <- as_factor(AERNASE.clin.hdac9$MedHx_CVD)
AERNASE.clin.hdac9$DiabetesStatus <- as_factor(AERNASE.clin.hdac9$DiabetesStatus)
study.samplesize = nrow(AERNASE.clin.hdac9) # study.samplesize is an expected variable in the GLM.BIN() GLM.CON() functions
TRAITS.TARGET.RANK = c("HDAC9")
TRAITS.CON.RANK = c("MAC_rankNorm", "SMC_rankNorm", "MAC_SMC_ratio_rank", "VesselDensity_rankNorm")
TRAITS.BIN = c("CalcificationPlaque", "CollagenPlaque", "Fat10Perc", "IPH",
"MAC_binned", "SMC_binned")
# "Hospital",
# "Age", "Gender",
# "TC_final", "LDL_final", "HDL_final", "TG_final",
# "systolic", "diastoli", "GFR_MDRD", "BMI",
# "KDOQI", "BMI_WHO",
# "SmokerCurrent", "eCigarettes", "ePackYearsSmoking",
# "DiabetesStatus", "Hypertension.composite",
# "Hypertension.drugs", "Med.anticoagulants", "Med.all.antiplatelet", "Med.Statin.LLD",
# "Stroke_Dx", "sympt", "Symptoms.5G", "restenos",
# "EP_composite", "EP_composite_time",
# "macmean0", "smcmean0", "Macrophages.bin", "SMC.bin",
# "neutrophils", "Mast_cells_plaque",
# "IPH.bin", "vessel_density_averaged",
# "Calc.bin", "Collagen.bin",
# "Fat.bin_10", "Fat.bin_40", "OverallPlaquePhenotype",
# "IL6_pg_ug_2015", "MCP1_pg_ug_2015",
# "QC2018_FILTER", "CHIP", "SAMPLE_TYPE",
# "CAD_history", "Stroke_history", "Peripheral.interv",
# "stenose"
# 1. Age (continuous in 1-year increment). [Age]
# 2. Sex (male vs. female). [Gender]
# 3. Presence of hypertension at baseline (defined either as history of hypertension, SBP ≥140 mm Hg, DBP ≥90 mm Hg, or prescription of antihypertensive medications). [Hypertension.composite]
# 4. Presence of diabetes mellitus at baseline (defined either as a history of diabetes, administration of glucose lowering medication, HbA1c ≥6.5%, fasting glucose ≥126 mg/dl, .or random glucose levels ≥200 mg/dl). [DiabetesStatus]
# 5. Smoking (current, ex-, never). [SmokerCurrent]
# 6. LDL-C levels (continuous). [LDL_final]
# 7. Use of lipid-lowering drugs. [Med.Statin.LLD]
# 8. Use of antiplatelet drugs. [Med.all.antiplatelet]
# 9. eGFR (continuous). [GFR_MDRD]
# 10. BMI (continuous). [BMI]
# 11. History of cardiovascular disease (stroke, coronary artery disease, peripheral artery disease). [MedHx_CVD] combinatino of: [CAD_history, Stroke_history, Peripheral.interv]
# 12. Level of stenosis (50-70% vs. 70-99%). [stenose]
# Models
# Model 1: adjusted for age and sex
# Model 2: adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis,
AERNASE.clin.hdac9$ORdate_epoch <- as.numeric(AERNASE.clin.hdac9$dateok)
AERNASE.clin.hdac9$ORdate_year <- AERNASE.clin.hdac9$ORyear
cat("Summary of 'year of surgery' in 'epoch' (); coded as `numeric()`\n")Summary of 'year of surgery' in 'epoch' (); coded as `numeric()`summary(AERNASE.clin.hdac9$ORdate_epoch) Min. 1st Qu. Median Mean 3rd Qu. Max.
11770 12759 13440 13553 14287 16807 cat("\nSummary of 'year of surgery' in 'years' (); coded as `factor()`\n")
Summary of 'year of surgery' in 'years' (); coded as `factor()`table(AERNASE.clin.hdac9$ORdate_year)
No data available/missing 2002 2003 2004 2005 2006
0 32 62 66 82 85
2007 2008 2009 2010 2011 2012
67 63 67 34 31 22
2013 2014 2015 2016 2017 2018
5 3 3 1 0 0
2019
0 # COVARIATES_M1 = c("Age", "Gender", "ORdate_year")
COVARIATES_M1 = c("Age", "Gender", "ORdate_epoch")
COVARIATES_M2 = c(COVARIATES_M1,
"Hypertension.composite", "DiabetesStatus",
"SmokerStatus",
# "SmokerCurrent",
"Med.Statin.LLD", "Med.all.antiplatelet",
"GFR_MDRD", "BMI",
# "CAD_history", "Stroke_history", "Peripheral.interv",
"MedHx_CVD",
"stenose")
str(AERNASE.clin.hdac9)'data.frame': 623 obs. of 74 variables:
$ STUDY_NUMBER : chr "ae1" "ae1026" "ae1029" "ae1032" ...
$ SampleType : chr "plaque" "plaque" "plaque" "plaque" ...
$ RNAseqType : chr "3' RNAseq" "3' RNAseq" "3' RNAseq" "3' RNAseq" ...
$ Hospital : Factor w/ 2 levels "St. Antonius, Nieuwegein",..: 1 2 2 2 2 2 2 1 2 2 ...
$ ORyear : Factor w/ 19 levels "No data available/missing",..: 2 5 5 5 5 5 6 3 5 7 ...
$ Artery_summary : Factor w/ 9 levels "No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA",..: 2 2 2 2 2 2 2 2 2 2 ...
$ Age : num 63 80 57 49 75 71 77 78 60 59 ...
$ Gender : Factor w/ 2 levels "female","male": 2 2 2 1 1 1 2 2 2 1 ...
$ TC_final : num 4.48 NA 4.2 NA 2.92 ...
$ LDL_final : num 2.94 NA NA NA 0.65 ...
$ HDL_final : num 1 NA NA NA 0.87 ...
$ TG_final : num 1.21 NA NA NA 3.07 ...
$ systolic : num 220 160 150 170 104 180 160 165 140 150 ...
$ diastoli : num 90 80 83 90 55 90 90 70 60 52 ...
$ GFR_MDRD : num 50.9 63.3 76.3 143.9 34.9 ...
$ BMI : num 28.6 26.5 21.6 22 30.8 ...
$ KDOQI : Factor w/ 6 levels "No data available/missing",..: 4 3 3 2 4 5 3 4 4 4 ...
$ BMI_WHO : Factor w/ 5 levels "No data available/missing",..: 4 4 3 3 5 4 4 3 4 3 ...
$ SmokerStatus : Factor w/ 3 levels "Current smoker",..: 2 NA 1 1 3 1 3 3 3 1 ...
$ AlcoholUse : Factor w/ 2 levels "No","Yes": 1 2 1 2 2 1 2 1 2 1 ...
$ Hypertension.selfreport : Factor w/ 3 levels "No data available/missing",..: 3 2 2 2 3 3 3 3 3 3 ...
$ Hypertension.selfreportdrug: Factor w/ 3 levels "No data available/missing",..: 3 2 2 2 3 2 3 3 3 3 ...
$ Hypertension.composite : Factor w/ 3 levels "No data available/missing",..: 3 3 2 2 3 3 3 3 3 3 ...
$ Hypertension.drugs : Factor w/ 3 levels "No data available/missing",..: 2 3 2 2 3 3 3 3 3 3 ...
$ Med.anticoagulants : Factor w/ 3 levels "No data available/missing",..: 2 3 2 2 2 3 2 3 2 2 ...
$ Med.all.antiplatelet : Factor w/ 3 levels "No data available/missing",..: 3 2 3 3 3 2 3 3 3 3 ...
$ Med.Statin.LLD : Factor w/ 3 levels "No data available/missing",..: 2 3 3 2 3 3 3 2 3 3 ...
$ Stroke_Dx : Factor w/ 3 levels "Missing","No stroke diagnosed",..: 2 2 2 2 2 2 2 2 2 2 ...
$ sympt : Factor w/ 18 levels "missing","Asymptomatic",..: 3 3 5 3 6 4 2 2 2 2 ...
$ Symptoms.5G : Factor w/ 6 levels "Asymptomatic",..: 6 6 5 6 2 5 1 1 1 1 ...
$ AsymptSympt : Factor w/ 3 levels "Asymptomatic",..: 3 3 3 3 2 3 1 1 1 1 ...
$ AsymptSympt2G : Factor w/ 2 levels "Asymptomatic",..: 2 2 2 2 2 2 1 1 1 1 ...
$ Symptoms.Update2G : Factor w/ 2 levels "Asymptomatic",..: 2 2 2 2 2 2 1 1 1 1 ...
$ Symptoms.Update3G : Factor w/ 3 levels "Asymptomatic",..: 2 2 2 2 2 2 1 1 1 1 ...
$ restenos : Factor w/ 4 levels "missing","de novo",..: 2 2 2 2 2 2 2 2 2 2 ...
$ stenose : Factor w/ 10 levels "missing","0-49%",..: 5 3 4 4 5 4 4 5 4 4 ...
$ CAD_history : Factor w/ 3 levels "Missing","No history CAD",..: 3 2 2 2 2 2 2 2 2 2 ...
$ PAOD : Factor w/ 3 levels "missing/no data",..: 3 3 2 2 2 2 2 2 2 3 ...
$ Peripheral.interv : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 2 ...
$ EP_composite : Factor w/ 3 levels "No data available.",..: 2 3 3 2 3 3 3 2 2 2 ...
$ EP_composite_time : num 12.258 0.208 0.685 1.288 2.603 ...
$ EP_major : Factor w/ 3 levels "No data available.",..: 2 2 2 2 3 3 3 2 2 2 ...
$ EP_major_time : num 12.258 0.647 2.989 1.288 2.603 ...
$ MAC_rankNorm : num -0.0695 -0.7504 0.9099 -0.3812 -1.0959 ...
$ SMC_rankNorm : num 1.526 0.58 0.674 1.26 0.341 ...
$ Macrophages.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 1 2 1 1 2 2 1 2 1 ...
$ SMC.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ Neutrophils_rankNorm : num NA 0.591 0.118 -1.618 -1.618 ...
$ MastCells_rankNorm : num NA 0.633 NA NA NA ...
$ IPH.bin : Factor w/ 2 levels "no","yes": 1 2 2 1 1 1 1 2 1 1 ...
$ VesselDensity_rankNorm : num -0.723 0.769 0.854 1.14 -0.861 ...
$ Calc.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 1 1 1 2 2 2 1 ...
$ Collagen.bin : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ Fat.bin_10 : Factor w/ 2 levels " <10%"," >10%": 2 2 2 2 1 1 2 1 2 1 ...
$ Fat.bin_40 : Factor w/ 2 levels "<40%",">40%": 1 1 1 1 1 1 2 1 1 1 ...
$ OverallPlaquePhenotype : Factor w/ 3 levels "atheromatous",..: 2 3 2 2 3 3 1 3 2 3 ...
$ Plaque_Vulnerability_Index : Factor w/ 6 levels "0","1","2","3",..: 3 3 4 2 1 2 5 2 3 1 ...
$ HDAC9 : int 24 11 7 5 13 46 14 0 30 12 ...
$ TWIST1 : int 1 13 5 0 1 1 3 1 19 2 ...
$ IL1B : int 2 4 2 1 1 2 5 2 43 2 ...
$ IL6 : int 0 1 0 0 0 0 0 0 5 0 ...
$ MAC_SMC_ratio : num -0.0455 -1.2938 1.3509 -0.3026 -3.2115 ...
$ MAC_SMC_ratio_rank : num -0.215 -0.888 0.722 -0.366 -1.331 ...
$ CalcificationPlaque : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 1 1 1 2 2 2 1 ...
$ CollagenPlaque : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ Fat10Perc : Factor w/ 2 levels " <10%"," >10%": 2 2 2 2 1 1 2 1 2 1 ...
$ MAC_binned : Factor w/ 2 levels "no/minor","moderate/heavy": 2 1 2 1 1 2 2 1 2 1 ...
$ SMC_binned : Factor w/ 2 levels "no/minor","moderate/heavy": 2 2 2 2 2 2 1 2 2 2 ...
$ IPH : Factor w/ 2 levels "no","yes": 1 2 2 1 1 1 1 2 1 1 ...
$ dateok : Date, format: "2002-03-26" "2005-09-08" "2005-09-15" "2005-09-29" ...
$ DiabetesStatus : Factor w/ 2 levels "Control (no Diabetes Dx/Med)",..: 1 1 1 1 2 1 1 1 1 1 ...
$ MedHx_CVD : Factor w/ 2 levels "yes","No": 1 2 1 2 2 1 2 2 2 1 ...
$ ORdate_epoch : num 11772 13034 13041 13055 13083 ...
$ ORdate_year : Factor w/ 19 levels "No data available/missing",..: 2 5 5 5 5 5 6 3 5 7 ...# COVARIATES_M3 = c(COVARIATES_M2, "LDL_final")
# COVARIATES_M4 = c(COVARIATES_M2, "hsCRP_plasma")In this model we correct for Age, Gender, and year of surgery.
Here we use the inverse-rank normalized data - visually this is more normally distributed.
detach("package:EnsDb.Hsapiens.v86", unload = TRUE)
detach("package:ensembldb", unload = TRUE)Analysis of continuous/quantitative plaque traits as a function of plaque HDAC9 expression levels.
library(dplyr)
library(MASS)
GLM.results <- data.frame(matrix(NA, ncol = 15, nrow = 0))
cat("Running linear regression...\n")Running linear regression...for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.CON.RANK)) {
TRAIT = TRAITS.CON.RANK[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# # for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
### univariate
# fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_year, data = currentDF)
fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_epoch, data = currentDF)
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 15, nrow = 0))
GLM.results.TEMP[1,] = GLM.CON(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing MAC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ 1, data = currentDF)
Coefficients:
(Intercept)
28.34
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-29.12 -17.76 -7.33 8.41 420.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.448488 20.604160 0.362 0.718
currentDF[, TRAIT] 0.406087 1.390528 0.292 0.770
Age 0.104075 0.148668 0.700 0.484
Gendermale -1.850604 3.049440 -0.607 0.544
ORdate_epoch 0.001114 0.001383 0.805 0.421
Residual standard error: 32.38 on 603 degrees of freedom
Multiple R-squared: 0.002608, Adjusted R-squared: -0.004008
F-statistic: 0.3942 on 4 and 603 DF, p-value: 0.8129
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_rankNorm
Effect size...............: 0.406087
Standard error............: 1.390528
Odds ratio (effect size)..: 1.501
Lower 95% CI..............: 0.098
Upper 95% CI..............: 22.909
T-value...................: 0.292038
P-value...................: 0.7703582
R^2.......................: 0.002608
Adjusted r^2..............: -0.004008
Sample size of AE DB......: 623
Sample size of model......: 608
Missing data %............: 2.407705
- processing SMC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT], data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT]
28.15 -4.19
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-33.64 -17.22 -7.40 7.66 420.33
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.3246475 20.8289123 1.024 0.30634
currentDF[, TRAIT] -4.1954910 1.4629836 -2.868 0.00428 **
Age 0.0282566 0.1497052 0.189 0.85035
Gendermale -2.8511622 3.0506642 -0.935 0.35037
ORdate_epoch 0.0005211 0.0013854 0.376 0.70696
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.19 on 601 degrees of freedom
Multiple R-squared: 0.01613, Adjusted R-squared: 0.009586
F-statistic: 2.464 on 4 and 601 DF, p-value: 0.04408
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_rankNorm
Effect size...............: -4.195491
Standard error............: 1.462984
Odds ratio (effect size)..: 0.015
Lower 95% CI..............: 0.001
Upper 95% CI..............: 0.265
T-value...................: -2.867764
P-value...................: 0.004278845
R^2.......................: 0.016134
Adjusted r^2..............: 0.009586
Sample size of AE DB......: 623
Sample size of model......: 606
Missing data %............: 2.728732
- processing MAC_SMC_ratio_rank
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT], data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT]
28.284 -2.673
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-33.52 -17.33 -7.38 7.23 418.71
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.347972 20.333310 0.263 0.7926
currentDF[, TRAIT] -2.823637 1.320233 -2.139 0.0329 *
Age 0.092945 0.148074 0.628 0.5304
Gendermale -2.417155 3.051326 -0.792 0.4286
ORdate_epoch 0.001362 0.001376 0.990 0.3226
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.29 on 601 degrees of freedom
Multiple R-squared: 0.0102, Adjusted R-squared: 0.003617
F-statistic: 1.549 on 4 and 601 DF, p-value: 0.1865
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_SMC_ratio_rank ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_SMC_ratio_rank
Effect size...............: -2.823637
Standard error............: 1.320233
Odds ratio (effect size)..: 0.059
Lower 95% CI..............: 0.004
Upper 95% CI..............: 0.79
T-value...................: -2.138741
P-value...................: 0.03285904
R^2.......................: 0.010204
Adjusted r^2..............: 0.003617
Sample size of AE DB......: 623
Sample size of model......: 606
Missing data %............: 2.728732
- processing VesselDensity_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ 1, data = currentDF)
Coefficients:
(Intercept)
28.7
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-29.91 -18.26 -7.72 8.32 420.46
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.289183 21.803793 0.243 0.808
currentDF[, TRAIT] -0.589271 1.475546 -0.399 0.690
Age 0.104402 0.156655 0.666 0.505
Gendermale -2.038897 3.207386 -0.636 0.525
ORdate_epoch 0.001328 0.001464 0.908 0.365
Residual standard error: 33.05 on 565 degrees of freedom
Multiple R-squared: 0.003443, Adjusted R-squared: -0.003612
F-statistic: 0.488 on 4 and 565 DF, p-value: 0.7446
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' VesselDensity_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: VesselDensity_rankNorm
Effect size...............: -0.589271
Standard error............: 1.475546
Odds ratio (effect size)..: 0.555
Lower 95% CI..............: 0.031
Upper 95% CI..............: 10.002
T-value...................: -0.399358
P-value...................: 0.6897803
R^2.......................: 0.003443
Adjusted r^2..............: -0.003612
Sample size of AE DB......: 623
Sample size of model......: 570
Missing data %............: 8.507223 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"T-value", "P-value", "r^2", "r^2_adj", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`T-value` <- as.numeric(GLM.results$`T-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2` <- as.numeric(GLM.results$`r^2`)
GLM.results$`r^2_adj` <- as.numeric(GLM.results$`r^2_adj`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
library(openxlsx)
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Con.Uni.",TRAIT_OF_INTEREST,".PlaquePhenotypes.RANK.MODEL1.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Con.Uni.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)Analysis of binary plaque traits as a function of plaque HDAC9 expression levels.
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.BIN)) {
TRAIT = TRAITS.BIN[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing CalcificationPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch
8.2993047 0.0058323 -0.3976708 -0.0005923
Degrees of Freedom: 614 Total (i.e. Null); 611 Residual
Null Deviance: 849.3
Residual Deviance: 793.8 AIC: 801.8
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.888 -1.118 0.774 1.044 1.708
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.973e+00 1.342e+00 5.939 2.87e-09 ***
currentDF[, TARGET] 5.798e-03 3.150e-03 1.840 0.0657 .
Age 5.704e-03 9.538e-03 0.598 0.5498
Gendermale -3.977e-01 1.993e-01 -1.995 0.0460 *
ORdate_epoch -5.970e-04 9.076e-05 -6.577 4.79e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 849.28 on 614 degrees of freedom
Residual deviance: 793.43 on 610 degrees of freedom
AIC: 803.43
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CalcificationPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CalcificationPlaque
Effect size...............: 0.005798
Standard error............: 0.00315
Odds ratio (effect size)..: 1.006
Lower 95% CI..............: 1
Upper 95% CI..............: 1.012
Z-value...................: 1.840355
P-value...................: 0.06571617
Hosmer and Lemeshow r^2...: 0.065755
Cox and Snell r^2.........: 0.086802
Nagelkerke's pseudo r^2...: 0.115944
Sample size of AE DB......: 623
Sample size of model......: 615
Missing data %............: 1.284109
- processing CollagenPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ 1, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept)
1.436
Degrees of Freedom: 613 Total (i.e. Null); 613 Residual
Null Deviance: 600.9
Residual Deviance: 600.9 AIC: 602.9
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9565 0.5931 0.6317 0.6687 1.0963
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.6386980 1.5759345 1.674 0.0941 .
currentDF[, TARGET] -0.0031744 0.0028088 -1.130 0.2584
Age 0.0030889 0.0115622 0.267 0.7893
Gendermale 0.1412298 0.2350181 0.601 0.5479
ORdate_epoch -0.0001054 0.0001045 -1.009 0.3131
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 600.95 on 613 degrees of freedom
Residual deviance: 598.35 on 609 degrees of freedom
AIC: 608.35
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CollagenPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CollagenPlaque
Effect size...............: -0.003174
Standard error............: 0.002809
Odds ratio (effect size)..: 0.997
Lower 95% CI..............: 0.991
Upper 95% CI..............: 1.002
Z-value...................: -1.130151
P-value...................: 0.2584124
Hosmer and Lemeshow r^2...: 0.004323
Cox and Snell r^2.........: 0.004222
Nagelkerke's pseudo r^2...: 0.006764
Sample size of AE DB......: 623
Sample size of model......: 614
Missing data %............: 1.444623
- processing Fat10Perc
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch
7.9906422 0.0194570 0.7733690 -0.0005761
Degrees of Freedom: 614 Total (i.e. Null); 611 Residual
Null Deviance: 671.8
Residual Deviance: 617.5 AIC: 625.5
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5274 0.2789 0.5799 0.7497 1.3251
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.2554211 1.5353858 4.725 2.30e-06 ***
currentDF[, TARGET] 0.0193811 0.0054251 3.573 0.000354 ***
Age 0.0127860 0.0110204 1.160 0.245962
Gendermale 0.7773361 0.2217286 3.506 0.000455 ***
ORdate_epoch -0.0005865 0.0001039 -5.646 1.65e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 671.77 on 614 degrees of freedom
Residual deviance: 616.15 on 610 degrees of freedom
AIC: 626.15
Number of Fisher Scoring iterations: 5
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' Fat10Perc ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: Fat10Perc
Effect size...............: 0.019381
Standard error............: 0.005425
Odds ratio (effect size)..: 1.02
Lower 95% CI..............: 1.009
Upper 95% CI..............: 1.03
Z-value...................: 3.572511
P-value...................: 0.0003535745
Hosmer and Lemeshow r^2...: 0.082805
Cox and Snell r^2.........: 0.08648
Nagelkerke's pseudo r^2...: 0.13013
Sample size of AE DB......: 623
Sample size of model......: 615
Missing data %............: 1.284109
- processing IPH
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + ORdate_epoch,
family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) Gendermale ORdate_epoch
7.8612128 0.5675364 -0.0005752
Degrees of Freedom: 613 Total (i.e. Null); 611 Residual
Null Deviance: 817.1
Residual Deviance: 768.2 AIC: 774.2
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9080 -1.1993 0.7252 0.9502 1.6768
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.555e+00 1.351e+00 5.591 2.26e-08 ***
currentDF[, TARGET] 1.453e-03 2.864e-03 0.507 0.61187
Age 4.985e-03 9.737e-03 0.512 0.60868
Gendermale 5.710e-01 1.982e-01 2.881 0.00396 **
ORdate_epoch -5.810e-04 9.146e-05 -6.353 2.12e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 817.10 on 613 degrees of freedom
Residual deviance: 767.62 on 609 degrees of freedom
AIC: 777.62
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' IPH ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: IPH
Effect size...............: 0.001453
Standard error............: 0.002864
Odds ratio (effect size)..: 1.001
Lower 95% CI..............: 0.996
Upper 95% CI..............: 1.007
Z-value...................: 0.507409
P-value...................: 0.6118676
Hosmer and Lemeshow r^2...: 0.060548
Cox and Snell r^2.........: 0.077415
Nagelkerke's pseudo r^2...: 0.105222
Sample size of AE DB......: 623
Sample size of model......: 614
Missing data %............: 1.444623
- processing MAC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale
-0.175134 0.004114 0.439912
Degrees of Freedom: 611 Total (i.e. Null); 609 Residual
Null Deviance: 837.4
Residual Deviance: 830 AIC: 836
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7896 -1.2867 0.9744 1.0453 1.2921
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.397e-01 1.255e+00 0.669 0.5033
currentDF[, TARGET] 4.233e-03 2.971e-03 1.425 0.1542
Age 1.051e-04 9.257e-03 0.011 0.9909
Gendermale 4.508e-01 1.898e-01 2.374 0.0176 *
ORdate_epoch -7.641e-05 8.404e-05 -0.909 0.3632
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 837.39 on 611 degrees of freedom
Residual deviance: 829.18 on 607 degrees of freedom
AIC: 839.18
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_binned
Effect size...............: 0.004233
Standard error............: 0.002971
Odds ratio (effect size)..: 1.004
Lower 95% CI..............: 0.998
Upper 95% CI..............: 1.01
Z-value...................: 1.424928
P-value...................: 0.154178
Hosmer and Lemeshow r^2...: 0.009808
Cox and Snell r^2.........: 0.013331
Nagelkerke's pseudo r^2...: 0.017883
Sample size of AE DB......: 623
Sample size of model......: 612
Missing data %............: 1.76565
- processing SMC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Age
2.436868 -0.004074 -0.022542
Degrees of Freedom: 611 Total (i.e. Null); 609 Residual
Null Deviance: 764.5
Residual Deviance: 756.7 AIC: 762.7
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7952 -1.4102 0.8040 0.8926 1.6860
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.992e+00 1.346e+00 2.222 0.0263 *
currentDF[, TARGET] -4.145e-03 2.702e-03 -1.534 0.1249
Age -2.222e-02 1.007e-02 -2.207 0.0273 *
Gendermale -2.770e-01 2.097e-01 -1.321 0.1866
ORdate_epoch -2.685e-05 8.954e-05 -0.300 0.7643
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 764.49 on 611 degrees of freedom
Residual deviance: 754.75 on 607 degrees of freedom
AIC: 764.75
Number of Fisher Scoring iterations: 4
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_binned
Effect size...............: -0.004145
Standard error............: 0.002702
Odds ratio (effect size)..: 0.996
Lower 95% CI..............: 0.991
Upper 95% CI..............: 1.001
Z-value...................: -1.534329
P-value...................: 0.1249488
Hosmer and Lemeshow r^2...: 0.012743
Cox and Snell r^2.........: 0.015792
Nagelkerke's pseudo r^2...: 0.022141
Sample size of AE DB......: 623
Sample size of model......: 612
Missing data %............: 1.76565 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Uni.",TRAIT_OF_INTEREST,".PlaquePhenotypes.RANK.MODEL1.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Bin.Uni.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)In this model we correct for Age, Gender, year of surgery, Hypertension status, Diabetes status, current smoker status, lipid-lowering drugs (LLDs), antiplatelet medication, eGFR (MDRD), BMI, MedHx_CVD (combination of CAD history, stroke history, and peripheral interventions), and stenosis.
Here we use the inverse-rank normalized data - visually
this is more normally distributed.
Analysis of continuous/quantitative plaque traits as a function of plaque HDAC9 expression levels.
GLM.results <- data.frame(matrix(NA, ncol = 15, nrow = 0))
cat("Running linear regression...\n")Running linear regression...for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.CON.RANK)) {
TRAIT = TRAITS.CON.RANK[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
### univariate
# fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF)
fit <- lm(currentDF[,TARGET] ~ currentDF[,TRAIT] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF)
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 15, nrow = 0))
GLM.results.TEMP[1,] = GLM.CON(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing MAC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ 1, data = currentDF)
Coefficients:
(Intercept)
28.59
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-32.79 -18.09 -6.73 7.70 414.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.091240 38.249933 -0.264 0.792
currentDF[, TRAIT] 0.011511 1.590977 0.007 0.994
Age 0.161680 0.195565 0.827 0.409
Gendermale 0.097035 3.568556 0.027 0.978
ORdate_epoch 0.001094 0.001635 0.669 0.504
Hypertension.compositeyes 2.904916 4.567294 0.636 0.525
DiabetesStatusDiabetes -5.303613 3.776907 -1.404 0.161
SmokerStatusEx-smoker -0.584559 3.421371 -0.171 0.864
SmokerStatusNever smoked -3.216563 4.525757 -0.711 0.478
Med.Statin.LLDyes -0.247286 3.648425 -0.068 0.946
Med.all.antiplateletyes 4.630310 5.206895 0.889 0.374
GFR_MDRD -0.021379 0.078478 -0.272 0.785
BMI -0.007456 0.407308 -0.018 0.985
MedHx_CVDNo 0.974358 3.120792 0.312 0.755
stenose50-70% 12.645258 24.734342 0.511 0.609
stenose70-90% 7.002051 24.074940 0.291 0.771
stenose90-99% 11.622272 24.096124 0.482 0.630
stenose100% (Occlusion) 7.828567 29.539505 0.265 0.791
Residual standard error: 33.64 on 507 degrees of freedom
Multiple R-squared: 0.01534, Adjusted R-squared: -0.01768
F-statistic: 0.4646 on 17 and 507 DF, p-value: 0.9676
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_rankNorm
Effect size...............: 0.011511
Standard error............: 1.590977
Odds ratio (effect size)..: 1.012
Lower 95% CI..............: 0.045
Upper 95% CI..............: 22.87
T-value...................: 0.007235
P-value...................: 0.9942299
R^2.......................: 0.015338
Adjusted r^2..............: -0.017678
Sample size of AE DB......: 623
Sample size of model......: 525
Missing data %............: 15.73034
- processing SMC_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT], data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT]
28.427 -4.472
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-39.66 -17.70 -6.56 6.78 413.59
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.7847801 38.2114545 -0.021 0.98362
currentDF[, TRAIT] -4.3670090 1.6496545 -2.647 0.00837 **
Age 0.1153929 0.1952382 0.591 0.55476
Gendermale -1.4107994 3.5852149 -0.394 0.69411
ORdate_epoch 0.0005486 0.0016237 0.338 0.73561
Hypertension.compositeyes 2.5480563 4.5364002 0.562 0.57458
DiabetesStatusDiabetes -4.4946625 3.7607644 -1.195 0.23259
SmokerStatusEx-smoker -0.4860758 3.4024177 -0.143 0.88646
SmokerStatusNever smoked -3.7090935 4.4983490 -0.825 0.41002
Med.Statin.LLDyes -0.5002184 3.6311515 -0.138 0.89049
Med.all.antiplateletyes 4.1991048 5.1746758 0.811 0.41748
GFR_MDRD -0.0021769 0.0781408 -0.028 0.97779
BMI -0.0058534 0.4045985 -0.014 0.98846
MedHx_CVDNo 1.5037054 3.1064665 0.484 0.62856
stenose50-70% 13.3478962 24.5399585 0.544 0.58673
stenose70-90% 8.2364390 23.8706266 0.345 0.73020
stenose90-99% 13.1154314 23.9133841 0.548 0.58362
stenose100% (Occlusion) 9.3401996 29.3473745 0.318 0.75042
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33.41 on 506 degrees of freedom
Multiple R-squared: 0.02898, Adjusted R-squared: -0.003638
F-statistic: 0.8885 on 17 and 506 DF, p-value: 0.5881
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_rankNorm
Effect size...............: -4.367009
Standard error............: 1.649655
Odds ratio (effect size)..: 0.013
Lower 95% CI..............: 0.001
Upper 95% CI..............: 0.322
T-value...................: -2.647226
P-value...................: 0.008368544
R^2.......................: 0.028985
Adjusted r^2..............: -0.003638
Sample size of AE DB......: 623
Sample size of model......: 524
Missing data %............: 15.89085
- processing MAC_SMC_ratio_rank
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + DiabetesStatus,
data = currentDF)
Coefficients:
(Intercept) currentDF[, TRAIT] DiabetesStatusDiabetes
29.591 -3.515 -5.360
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-36.46 -16.97 -6.66 8.36 410.03
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -20.439631 38.126367 -0.536 0.59212
currentDF[, TRAIT] -4.066568 1.522841 -2.670 0.00782 **
Age 0.164831 0.194296 0.848 0.39664
Gendermale -0.561417 3.549429 -0.158 0.87439
ORdate_epoch 0.001597 0.001614 0.989 0.32300
Hypertension.compositeyes 3.311944 4.537894 0.730 0.46582
DiabetesStatusDiabetes -5.674795 3.754089 -1.512 0.13125
SmokerStatusEx-smoker -0.964921 3.402408 -0.284 0.77683
SmokerStatusNever smoked -4.245343 4.510008 -0.941 0.34699
Med.Statin.LLDyes 0.004525 3.625219 0.001 0.99900
Med.all.antiplateletyes 5.475340 5.184680 1.056 0.29144
GFR_MDRD -0.020198 0.077888 -0.259 0.79549
BMI -0.069206 0.405156 -0.171 0.86444
MedHx_CVDNo 1.327313 3.103312 0.428 0.66905
stenose50-70% 16.043342 24.567862 0.653 0.51404
stenose70-90% 11.067584 23.912034 0.463 0.64367
stenose90-99% 16.311118 23.965602 0.681 0.49643
stenose100% (Occlusion) 11.942968 29.378022 0.407 0.68453
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33.41 on 506 degrees of freedom
Multiple R-squared: 0.02922, Adjusted R-squared: -0.003397
F-statistic: 0.8958 on 17 and 506 DF, p-value: 0.5792
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_SMC_ratio_rank ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_SMC_ratio_rank
Effect size...............: -4.066568
Standard error............: 1.522841
Odds ratio (effect size)..: 0.017
Lower 95% CI..............: 0.001
Upper 95% CI..............: 0.339
T-value...................: -2.670383
P-value...................: 0.007820257
R^2.......................: 0.029218
Adjusted r^2..............: -0.003397
Sample size of AE DB......: 623
Sample size of model......: 524
Missing data %............: 15.89085
- processing VesselDensity_rankNorm
Call:
lm(formula = currentDF[, TARGET] ~ DiabetesStatus, data = currentDF)
Coefficients:
(Intercept) DiabetesStatusDiabetes
30.094 -6.184
Call:
lm(formula = currentDF[, TARGET] ~ currentDF[, TRAIT] + Age +
Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = currentDF)
Residuals:
Min 1Q Median 3Q Max
-33.84 -18.19 -6.70 8.32 414.11
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -13.375880 40.431373 -0.331 0.741
currentDF[, TRAIT] -0.198401 1.717763 -0.115 0.908
Age 0.146617 0.207997 0.705 0.481
Gendermale 0.040816 3.759758 0.011 0.991
ORdate_epoch 0.001486 0.001739 0.855 0.393
Hypertension.compositeyes 0.910538 4.861182 0.187 0.851
DiabetesStatusDiabetes -6.394737 4.140531 -1.544 0.123
SmokerStatusEx-smoker -0.647605 3.642397 -0.178 0.859
SmokerStatusNever smoked -2.855860 4.856549 -0.588 0.557
Med.Statin.LLDyes -1.765814 3.849786 -0.459 0.647
Med.all.antiplateletyes 4.990184 5.525657 0.903 0.367
GFR_MDRD -0.038742 0.083117 -0.466 0.641
BMI 0.090860 0.435965 0.208 0.835
MedHx_CVDNo 1.807646 3.312055 0.546 0.585
stenose50-70% 12.237496 25.415228 0.482 0.630
stenose70-90% 7.626824 24.609474 0.310 0.757
stenose90-99% 12.192396 24.646142 0.495 0.621
stenose100% (Occlusion) 7.967067 30.242214 0.263 0.792
stenose50-99% 3.506790 42.857537 0.082 0.935
Residual standard error: 34.39 on 473 degrees of freedom
Multiple R-squared: 0.0167, Adjusted R-squared: -0.02072
F-statistic: 0.4463 on 18 and 473 DF, p-value: 0.9771
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' VesselDensity_rankNorm ' .
Collecting data.
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: VesselDensity_rankNorm
Effect size...............: -0.198401
Standard error............: 1.717763
Odds ratio (effect size)..: 0.82
Lower 95% CI..............: 0.028
Upper 95% CI..............: 23.77
T-value...................: -0.1155
P-value...................: 0.9080981
R^2.......................: 0.0167
Adjusted r^2..............: -0.020719
Sample size of AE DB......: 623
Sample size of model......: 492
Missing data %............: 21.02729 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"T-value", "P-value", "r^2", "r^2_adj", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`T-value` <- as.numeric(GLM.results$`T-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2` <- as.numeric(GLM.results$`r^2`)
GLM.results$`r^2_adj` <- as.numeric(GLM.results$`r^2_adj`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
library(openxlsx)
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Con.Multi.",TRAIT_OF_INTEREST,".PlaquePhenotypes.RANK.MODEL2.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Con.Multi.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)Analysis of binary plaque traits as a function of plaque MCP1 levels.
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
for (trait in 1:length(TRAITS.BIN)) {
TRAIT = TRAITS.BIN[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
Analysis of HDAC9.
- processing CalcificationPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch + Hypertension.composite, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch Hypertension.compositeyes
6.8650431 0.0056556 -0.3311542 -0.0005261 0.4837455
Degrees of Freedom: 529 Total (i.e. Null); 525 Residual
Null Deviance: 733.8
Residual Deviance: 691.8 AIC: 701.8
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7749 -1.0842 0.7152 1.0645 1.7710
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.263e+00 2.394e+00 3.034 0.00241 **
currentDF[, TARGET] 5.323e-03 3.394e-03 1.568 0.11686
Age -4.408e-03 1.212e-02 -0.364 0.71604
Gendermale -3.524e-01 2.227e-01 -1.583 0.11346
ORdate_epoch -4.677e-04 1.007e-04 -4.645 3.41e-06 ***
Hypertension.compositeyes 4.955e-01 2.852e-01 1.737 0.08231 .
DiabetesStatusDiabetes -8.297e-02 2.322e-01 -0.357 0.72089
SmokerStatusEx-smoker -2.505e-02 2.124e-01 -0.118 0.90611
SmokerStatusNever smoked -2.432e-01 2.810e-01 -0.865 0.38683
Med.Statin.LLDyes -1.093e-01 2.254e-01 -0.485 0.62760
Med.all.antiplateletyes -4.505e-02 3.217e-01 -0.140 0.88862
GFR_MDRD -4.935e-03 4.848e-03 -1.018 0.30866
BMI -4.248e-04 2.524e-02 -0.017 0.98657
MedHx_CVDNo 2.581e-02 1.940e-01 0.133 0.89418
stenose50-70% -4.494e-01 1.542e+00 -0.291 0.77079
stenose70-90% -5.059e-01 1.499e+00 -0.338 0.73574
stenose90-99% -1.327e-01 1.502e+00 -0.088 0.92956
stenose100% (Occlusion) 6.675e-01 1.907e+00 0.350 0.72638
stenose50-99% -1.218e+01 5.354e+02 -0.023 0.98185
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 733.82 on 529 degrees of freedom
Residual deviance: 684.37 on 511 degrees of freedom
AIC: 722.37
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CalcificationPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CalcificationPlaque
Effect size...............: 0.005323
Standard error............: 0.003394
Odds ratio (effect size)..: 1.005
Lower 95% CI..............: 0.999
Upper 95% CI..............: 1.012
Z-value...................: 1.568105
P-value...................: 0.1168566
Hosmer and Lemeshow r^2...: 0.067384
Cox and Snell r^2.........: 0.089078
Nagelkerke's pseudo r^2...: 0.118839
Sample size of AE DB......: 623
Sample size of model......: 530
Missing data %............: 14.92777
- processing CollagenPlaque
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ SmokerStatus +
Med.all.antiplatelet, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) SmokerStatusEx-smoker SmokerStatusNever smoked Med.all.antiplateletyes
1.08806 -0.08434 -0.66574 0.61326
Degrees of Freedom: 528 Total (i.e. Null); 525 Residual
Null Deviance: 510
Residual Deviance: 501.9 AIC: 509.9
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1896 0.4820 0.5659 0.6630 1.1019
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.320e+01 1.009e+03 0.013 0.9896
currentDF[, TARGET] -3.871e-03 3.017e-03 -1.283 0.1994
Age 1.035e-02 1.506e-02 0.687 0.4919
Gendermale -1.286e-02 2.757e-01 -0.047 0.9628
ORdate_epoch 8.164e-05 1.247e-04 0.655 0.5126
Hypertension.compositeyes -1.573e-01 3.621e-01 -0.434 0.6639
DiabetesStatusDiabetes 1.813e-04 2.959e-01 0.001 0.9995
SmokerStatusEx-smoker -1.208e-01 2.755e-01 -0.438 0.6611
SmokerStatusNever smoked -7.115e-01 3.272e-01 -2.174 0.0297 *
Med.Statin.LLDyes -1.921e-01 2.866e-01 -0.670 0.5026
Med.all.antiplateletyes 7.816e-01 3.579e-01 2.184 0.0290 *
GFR_MDRD -1.859e-03 6.119e-03 -0.304 0.7613
BMI 2.833e-02 3.252e-02 0.871 0.3836
MedHx_CVDNo -1.020e-01 2.415e-01 -0.422 0.6728
stenose50-70% -1.506e+01 1.009e+03 -0.015 0.9881
stenose70-90% -1.429e+01 1.009e+03 -0.014 0.9887
stenose90-99% -1.403e+01 1.009e+03 -0.014 0.9889
stenose100% (Occlusion) 2.998e-01 1.233e+03 0.000 0.9998
stenose50-99% -6.314e-01 1.771e+03 0.000 0.9997
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 510.02 on 528 degrees of freedom
Residual deviance: 490.30 on 510 degrees of freedom
AIC: 528.3
Number of Fisher Scoring iterations: 14
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' CollagenPlaque ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: CollagenPlaque
Effect size...............: -0.003871
Standard error............: 0.003017
Odds ratio (effect size)..: 0.996
Lower 95% CI..............: 0.99
Upper 95% CI..............: 1.002
Z-value...................: -1.283273
P-value...................: 0.1993966
Hosmer and Lemeshow r^2...: 0.038651
Cox and Snell r^2.........: 0.036579
Nagelkerke's pseudo r^2...: 0.059124
Sample size of AE DB......: 623
Sample size of model......: 529
Missing data %............: 15.08828
- processing Fat10Perc
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Gendermale ORdate_epoch
8.1301039 0.0203375 0.6845070 -0.0005857
Degrees of Freedom: 529 Total (i.e. Null); 526 Residual
Null Deviance: 583.7
Residual Deviance: 535.6 AIC: 543.6
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.63164 0.09199 0.55976 0.74799 1.36449
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.200e+01 6.125e+02 0.036 0.97135
currentDF[, TARGET] 1.990e-02 6.053e-03 3.288 0.00101 **
Age 2.015e-02 1.432e-02 1.407 0.15930
Gendermale 7.426e-01 2.527e-01 2.939 0.00329 **
ORdate_epoch -6.930e-04 1.205e-04 -5.750 8.93e-09 ***
Hypertension.compositeyes -1.126e-02 3.375e-01 -0.033 0.97339
DiabetesStatusDiabetes -1.478e-01 2.727e-01 -0.542 0.58786
SmokerStatusEx-smoker -1.207e-01 2.530e-01 -0.477 0.63329
SmokerStatusNever smoked -1.047e-01 3.328e-01 -0.315 0.75298
Med.Statin.LLDyes 1.378e-01 2.720e-01 0.507 0.61249
Med.all.antiplateletyes 4.566e-01 3.636e-01 1.256 0.20928
GFR_MDRD 4.535e-03 5.812e-03 0.780 0.43524
BMI -3.618e-02 2.910e-02 -1.243 0.21382
MedHx_CVDNo 6.926e-02 2.296e-01 0.302 0.76288
stenose50-70% -1.266e+01 6.125e+02 -0.021 0.98351
stenose70-90% -1.353e+01 6.125e+02 -0.022 0.98237
stenose90-99% -1.378e+01 6.125e+02 -0.022 0.98205
stenose100% (Occlusion) -1.429e+01 6.125e+02 -0.023 0.98138
stenose50-99% -2.813e+01 1.074e+03 -0.026 0.97911
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 583.68 on 529 degrees of freedom
Residual deviance: 521.56 on 511 degrees of freedom
AIC: 559.56
Number of Fisher Scoring iterations: 13
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' Fat10Perc ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: Fat10Perc
Effect size...............: 0.019905
Standard error............: 0.006053
Odds ratio (effect size)..: 1.02
Lower 95% CI..............: 1.008
Upper 95% CI..............: 1.032
Z-value...................: 3.28842
P-value...................: 0.001007514
Hosmer and Lemeshow r^2...: 0.10643
Cox and Snell r^2.........: 0.110601
Nagelkerke's pseudo r^2...: 0.165681
Sample size of AE DB......: 623
Sample size of model......: 530
Missing data %............: 14.92777
- processing IPH
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + ORdate_epoch +
MedHx_CVD, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) Gendermale ORdate_epoch MedHx_CVDNo
7.6697898 0.3451109 -0.0005425 -0.2775456
Degrees of Freedom: 528 Total (i.e. Null); 525 Residual
Null Deviance: 705.4
Residual Deviance: 667.2 AIC: 675.2
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0497 -1.1745 0.7199 0.9686 1.7180
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.208e+01 5.807e+02 0.038 0.9697
currentDF[, TARGET] 1.632e-03 2.994e-03 0.545 0.5857
Age 2.692e-03 1.244e-02 0.216 0.8288
Gendermale 3.744e-01 2.241e-01 1.670 0.0949 .
ORdate_epoch -5.413e-04 1.034e-04 -5.237 1.63e-07 ***
Hypertension.compositeyes -2.178e-01 2.913e-01 -0.748 0.4547
DiabetesStatusDiabetes -1.030e-01 2.386e-01 -0.432 0.6660
SmokerStatusEx-smoker 1.670e-02 2.178e-01 0.077 0.9389
SmokerStatusNever smoked 1.959e-01 2.888e-01 0.678 0.4977
Med.Statin.LLDyes -2.582e-01 2.360e-01 -1.094 0.2739
Med.all.antiplateletyes 2.576e-01 3.278e-01 0.786 0.4320
GFR_MDRD -2.910e-03 4.981e-03 -0.584 0.5591
BMI 4.573e-03 2.580e-02 0.177 0.8593
MedHx_CVDNo -3.563e-01 1.972e-01 -1.807 0.0707 .
stenose50-70% -1.444e+01 5.807e+02 -0.025 0.9802
stenose70-90% -1.448e+01 5.807e+02 -0.025 0.9801
stenose90-99% -1.436e+01 5.807e+02 -0.025 0.9803
stenose100% (Occlusion) -1.517e+01 5.807e+02 -0.026 0.9792
stenose50-99% -2.844e+01 1.057e+03 -0.027 0.9785
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 705.43 on 528 degrees of freedom
Residual deviance: 658.40 on 510 degrees of freedom
AIC: 696.4
Number of Fisher Scoring iterations: 13
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' IPH ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: IPH
Effect size...............: 0.001632
Standard error............: 0.002994
Odds ratio (effect size)..: 1.002
Lower 95% CI..............: 0.996
Upper 95% CI..............: 1.008
Z-value...................: 0.545074
P-value...................: 0.5857023
Hosmer and Lemeshow r^2...: 0.066668
Cox and Snell r^2.........: 0.085065
Nagelkerke's pseudo r^2...: 0.115507
Sample size of AE DB......: 623
Sample size of model......: 529
Missing data %............: 15.08828
- processing MAC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ Gender + BMI, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) Gendermale BMI
0.96605 0.42042 -0.03957
Degrees of Freedom: 526 Total (i.e. Null); 524 Residual
Null Deviance: 723.5
Residual Deviance: 716.7 AIC: 722.7
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7916 -1.2337 0.9035 1.0567 1.5594
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.901e+00 2.293e+00 0.829 0.4071
currentDF[, TARGET] 3.554e-03 3.066e-03 1.159 0.2463
Age -5.878e-04 1.186e-02 -0.050 0.9605
Gendermale 4.628e-01 2.158e-01 2.145 0.0320 *
ORdate_epoch -9.757e-05 9.660e-05 -1.010 0.3125
Hypertension.compositeyes 1.367e-01 2.746e-01 0.498 0.6188
DiabetesStatusDiabetes -1.371e-01 2.298e-01 -0.597 0.5508
SmokerStatusEx-smoker 9.321e-02 2.078e-01 0.449 0.6538
SmokerStatusNever smoked 1.328e-01 2.745e-01 0.484 0.6286
Med.Statin.LLDyes 2.739e-01 2.198e-01 1.246 0.2127
Med.all.antiplateletyes -1.075e-01 3.174e-01 -0.339 0.7348
GFR_MDRD -2.536e-03 4.770e-03 -0.532 0.5950
BMI -4.501e-02 2.487e-02 -1.809 0.0704 .
MedHx_CVDNo -5.758e-02 1.895e-01 -0.304 0.7612
stenose50-70% 4.294e-02 1.473e+00 0.029 0.9767
stenose70-90% 5.184e-01 1.431e+00 0.362 0.7171
stenose90-99% 3.018e-01 1.433e+00 0.211 0.8332
stenose100% (Occlusion) -1.287e-01 1.760e+00 -0.073 0.9417
stenose50-99% -1.303e+01 5.354e+02 -0.024 0.9806
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 723.50 on 526 degrees of freedom
Residual deviance: 706.23 on 508 degrees of freedom
AIC: 744.23
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' MAC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: MAC_binned
Effect size...............: 0.003554
Standard error............: 0.003066
Odds ratio (effect size)..: 1.004
Lower 95% CI..............: 0.998
Upper 95% CI..............: 1.01
Z-value...................: 1.15947
P-value...................: 0.2462649
Hosmer and Lemeshow r^2...: 0.023877
Cox and Snell r^2.........: 0.032249
Nagelkerke's pseudo r^2...: 0.043193
Sample size of AE DB......: 623
Sample size of model......: 527
Missing data %............: 15.40931
- processing SMC_binned
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + MedHx_CVD, family = binomial(link = "logit"),
data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Age Gendermale MedHx_CVDNo
2.11537 -0.00535 -0.01568 -0.34967 0.43278
Degrees of Freedom: 527 Total (i.e. Null); 523 Residual
Null Deviance: 659
Residual Deviance: 645.5 AIC: 655.5
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0716 -1.2985 0.7195 0.8876 1.7432
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.770e-01 2.415e+00 -0.156 0.8759
currentDF[, TARGET] -6.093e-03 2.918e-03 -2.088 0.0368 *
Age -1.910e-02 1.309e-02 -1.459 0.1446
Gendermale -3.831e-01 2.439e-01 -1.571 0.1162
ORdate_epoch 1.153e-04 1.064e-04 1.083 0.2786
Hypertension.compositeyes 1.180e-01 2.978e-01 0.396 0.6918
DiabetesStatusDiabetes -1.661e-01 2.449e-01 -0.678 0.4975
SmokerStatusEx-smoker -4.962e-03 2.262e-01 -0.022 0.9825
SmokerStatusNever smoked -4.407e-01 2.911e-01 -1.514 0.1301
Med.Statin.LLDyes 5.922e-02 2.375e-01 0.249 0.8031
Med.all.antiplateletyes -2.639e-01 3.497e-01 -0.755 0.4504
GFR_MDRD -3.020e-03 5.139e-03 -0.588 0.5568
BMI 3.035e-02 2.771e-02 1.095 0.2734
MedHx_CVDNo 4.809e-01 2.101e-01 2.289 0.0221 *
stenose50-70% 2.059e-01 1.475e+00 0.140 0.8890
stenose70-90% 6.309e-01 1.433e+00 0.440 0.6597
stenose90-99% 1.122e+00 1.438e+00 0.781 0.4350
stenose100% (Occlusion) -5.214e-02 1.773e+00 -0.029 0.9765
stenose50-99% 1.380e+01 5.354e+02 0.026 0.9794
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 658.99 on 527 degrees of freedom
Residual deviance: 629.87 on 509 degrees of freedom
AIC: 667.87
Number of Fisher Scoring iterations: 12
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' SMC_binned ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: SMC_binned
Effect size...............: -0.006093
Standard error............: 0.002918
Odds ratio (effect size)..: 0.994
Lower 95% CI..............: 0.988
Upper 95% CI..............: 1
Z-value...................: -2.087993
P-value...................: 0.03679844
Hosmer and Lemeshow r^2...: 0.044176
Cox and Snell r^2.........: 0.053642
Nagelkerke's pseudo r^2...: 0.075241
Sample size of AE DB......: 623
Sample size of model......: 528
Missing data %............: 15.2488 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Multi.",TRAIT_OF_INTEREST,".PlaquePhenotypes.RANK.MODEL2.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Bin.Multi.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)We will perform a cross-sectional analysis between plaque HDAC9 expression levels and the ‘clinical status’ of the plaque in terms of presence of patients’ symptoms (symptomatic vs. asymptomatic). The symptoms of interest are:
In this model we correct for Age, Gender, and year of surgery.
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
TRAIT = "AsymptSympt"
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# + Hypertension.composite + DiabetesStatus + SmokerCurrent +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# CAD_history + Stroke_history + Peripheral.interv + stenose
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
Analysis of HDAC9.
- processing AsymptSympt
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Age
-1.21372 0.01044 0.04260
Degrees of Freedom: 622 Total (i.e. Null); 620 Residual
Null Deviance: 477.7
Residual Deviance: 463.5 AIC: 469.5
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.5580 0.3828 0.4791 0.5704 0.9491
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.961e+00 1.820e+00 -1.078 0.28117
currentDF[, TARGET] 9.840e-03 5.996e-03 1.641 0.10078
Age 4.313e-02 1.363e-02 3.164 0.00155 **
Gendermale -4.060e-01 3.070e-01 -1.322 0.18606
ORdate_epoch 7.727e-05 1.186e-04 0.652 0.51458
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 477.66 on 622 degrees of freedom
Residual deviance: 461.27 on 618 degrees of freedom
AIC: 471.27
Number of Fisher Scoring iterations: 5
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' AsymptSympt ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: AsymptSympt
Effect size...............: 0.00984
Standard error............: 0.005996
Odds ratio (effect size)..: 1.01
Lower 95% CI..............: 0.998
Upper 95% CI..............: 1.022
Z-value...................: 1.64107
P-value...................: 0.1007828
Hosmer and Lemeshow r^2...: 0.03432
Cox and Snell r^2.........: 0.02597
Nagelkerke's pseudo r^2...: 0.0485
Sample size of AE DB......: 623
Sample size of model......: 623
Missing data %............: 0 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Uni.",TRAIT_OF_INTEREST,".RANK.Symptoms.MODEL1.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Bin.Uni.Symptoms")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)In this model we correct for Age, Gender, Hypertension status, Diabetes status, current smoker status, lipid-lowering drugs (LLDs), antiplatelet medication, eGFR (MDRD), BMI, MedHx_CVD (combination of CAD history, stroke history, and peripheral interventions), and stenosis..
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)) {
TARGET = TRAITS.TARGET.RANK[target_of_interest]
cat(paste0("\nAnalysis of ",TARGET,".\n"))
TRAIT = "AsymptSympt"
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., TARGET, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF, family = binomial(link = "logit"))
fit <- glm(as.factor(currentDF[,TRAIT]) ~ currentDF[,TARGET] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF, family = binomial(link = "logit"))
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 16, nrow = 0))
GLM.results.TEMP[1,] = GLM.BIN(fit, "AERNASE.clin.hdac9", TARGET, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
Analysis of HDAC9.
- processing AsymptSympt
Call: glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + BMI, family = binomial(link = "logit"), data = currentDF)
Coefficients:
(Intercept) currentDF[, TARGET] Age BMI
0.583817 0.008055 0.037885 -0.052261
Degrees of Freedom: 536 Total (i.e. Null); 533 Residual
Null Deviance: 408
Residual Deviance: 396.3 AIC: 404.3
Call:
glm(formula = as.factor(currentDF[, TRAIT]) ~ currentDF[, TARGET] +
Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, family = binomial(link = "logit"),
data = currentDF)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6448 0.3364 0.4649 0.5673 0.9144
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.467e+01 1.690e+03 0.009 0.993
currentDF[, TARGET] 8.221e-03 6.461e-03 1.272 0.203
Age 4.453e-02 1.773e-02 2.512 0.012 *
Gendermale -2.256e-01 3.327e-01 -0.678 0.498
ORdate_epoch 7.762e-05 1.357e-04 0.572 0.567
Hypertension.compositeyes -3.389e-01 4.693e-01 -0.722 0.470
DiabetesStatusDiabetes 3.098e-02 3.300e-01 0.094 0.925
SmokerStatusEx-smoker -3.760e-01 3.154e-01 -1.192 0.233
SmokerStatusNever smoked -5.411e-01 3.989e-01 -1.357 0.175
Med.Statin.LLDyes -2.513e-01 3.538e-01 -0.710 0.477
Med.all.antiplateletyes -2.599e-01 5.583e-01 -0.466 0.642
GFR_MDRD 1.138e-03 7.188e-03 0.158 0.874
BMI -4.275e-02 3.442e-02 -1.242 0.214
MedHx_CVDNo 7.272e-02 2.819e-01 0.258 0.796
stenose50-70% -1.333e+01 1.690e+03 -0.008 0.994
stenose70-90% -1.483e+01 1.690e+03 -0.009 0.993
stenose90-99% -1.484e+01 1.690e+03 -0.009 0.993
stenose100% (Occlusion) -3.106e-01 2.066e+03 0.000 1.000
stenose50-99% -7.321e-01 2.935e+03 0.000 1.000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 408.04 on 536 degrees of freedom
Residual deviance: 384.85 on 518 degrees of freedom
AIC: 422.85
Number of Fisher Scoring iterations: 15
Analyzing in dataset ' AERNASE.clin.hdac9 ' the association of ' HDAC9 ' with ' AsymptSympt ' ...
Collecting data...
We have collected the following and summarize it in an object:
Dataset...................: AERNASE.clin.hdac9
Score/Exposure/biomarker..: HDAC9
Trait/outcome.............: AsymptSympt
Effect size...............: 0.008221
Standard error............: 0.006461
Odds ratio (effect size)..: 1.008
Lower 95% CI..............: 0.996
Upper 95% CI..............: 1.021
Z-value...................: 1.272457
P-value...................: 0.203211
Hosmer and Lemeshow r^2...: 0.056844
Cox and Snell r^2.........: 0.042274
Nagelkerke's pseudo r^2...: 0.079423
Sample size of AE DB......: 623
Sample size of model......: 537
Missing data %............: 13.80417 cat("Edit the column names...\n")Edit the column names...colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"Z-value", "P-value", "r^2_l", "r^2_cs", "r^2_nagelkerke", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")Correct the variable types...GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`Z-value` <- as.numeric(GLM.results$`Z-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2_l` <- as.numeric(GLM.results$`r^2_l`)
GLM.results$`r^2_cs` <- as.numeric(GLM.results$`r^2_cs`)
GLM.results$`r^2_nagelkerke` <- as.numeric(GLM.results$`r^2_nagelkerke`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)
# Save the data
cat("Writing results to Excel-file...\n")Writing results to Excel-file...### Univariate
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Bin.Multi.",TRAIT_OF_INTEREST,".RANK.Symptoms.MODEL2.xlsx"),
rowNames = FALSE, colNmes = TRUE, sheetName = "Bin.Multi.Symptoms")
# Removing intermediates
cat("Removing intermediate files...\n")Removing intermediate files...rm(TRAIT, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)For the longitudinal analyses of plaque HDAC9 expression levels and secondary cardiovascular events over a three-year follow-up period.
The primary outcome is defined as “a composite of fatal or
non-fatal myocardial infarction, fatal or non-fatal stroke, ruptured
aortic aneurysm, fatal cardiac failure, coronary or peripheral
interventions, leg amputation due to vascular causes, and cardiovascular
death”, i.e. major adverse cardiovascular events (MACE). Variable:
epmajor.3years, these include: - myocardial infarction (MI)
- cerebral infarction (CVA/stroke) - cardiovascular death (exact cause
to be investigated) - cerebral bleeding (CVA/stroke) - fatal myocardial
infarction (MI) - fatal cerebral infarction - fatal cerebral bleeding -
sudden death - fatal heart failure - fatal aneurysm rupture - other
cardiovascular death..
The secondary outcomes will be
epstroke.3years, these include:
epcoronary.3years, these include:
epcvdeath.3years,
these include:
We will use 3-year follow-up, but we will also calculate 30 days and 90 days follow-up ‘time-to-event’ variables. On average there are 365.25 days in a year. We can calculate 30-days and 90-days follow-up time based on the three years follow-up.
cutt.off.30days = (1/365.25) * 30
cutt.off.90days = (1/365.25) * 90
# Fix maximum FU of 30 and 90 days
AEDB.CEA <- AEDB.CEA %>%
mutate(
FU.cutt.off.30days = ifelse(max.followup <= cutt.off.30days, max.followup, cutt.off.30days),
FU.cutt.off.90days = ifelse(max.followup <= cutt.off.90days, max.followup, cutt.off.90days)
)
AEDB.temp <- subset(AEDB.CEA, select = c("STUDY_NUMBER", "Age", "Gender", "Hospital", "Artery_summary",
"max.followup",
"FU.cutt.off.3years",
"FU.cutt.off.30days",
"FU.cutt.off.90days"))
require(labelled)
AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
rm(AEDB.temp)
AERNASE.clin.hdac9 <- merge(AERNASE.clin.hdac9,
subset(AEDB.CEA, select = c("STUDY_NUMBER", "max.followup",
"FU.cutt.off.3years", "FU.cutt.off.30days", "FU.cutt.off.90days",
"ep_major_t_3years", "ep_stroke_t_3years", "ep_coronary_t_3years",
"ep_cvdeath_t_3years",
"epmajor.3years", "epstroke.3years", "epcoronary.3years", "epcvdeath.3years"
)),
by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE)
AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9, select = c("STUDY_NUMBER", "Age", "Gender", "Hospital", "Artery_summary",
"max.followup",
"FU.cutt.off.3years",
"FU.cutt.off.30days",
"FU.cutt.off.90days"))
require(labelled)
AERNASE.clin.hdac9.temp$Gender <- to_factor(AERNASE.clin.hdac9.temp$Gender)
AERNASE.clin.hdac9.temp$Hospital <- to_factor(AERNASE.clin.hdac9.temp$Hospital)
AERNASE.clin.hdac9.temp$Artery_summary <- to_factor(AERNASE.clin.hdac9.temp$Artery_summary)
DT::datatable(AERNASE.clin.hdac9.temp[1:10,], caption = "Excerpt of the whole AERNASE.clin.hdac9.", rownames = FALSE)
rm(AERNASE.clin.hdac9.temp)Here we will calculate the new 30- and 90-days follow-up of the events and their event-times of interest:
epmajor.3years)epstroke.3years)epcoronary.3years)epcvdeath.3years)avg_days_in_year = 365.25
cutt.off.30days.scaled <- cutt.off.30days * 365.25
cutt.off.90days.scaled <- cutt.off.90days * 365.25
# Event times
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>%
mutate(
ep_major_t_30days = ifelse(ep_major_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_major_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_stroke_t_30days = ifelse(ep_stroke_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_stroke_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_coronary_t_30days = ifelse(ep_coronary_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_coronary_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_cvdeath_t_30days = ifelse(ep_cvdeath_t_3years * avg_days_in_year <= cutt.off.30days.scaled,
ep_cvdeath_t_3years * avg_days_in_year, cutt.off.30days.scaled),
ep_major_t_90days = ifelse(ep_major_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_major_t_3years * avg_days_in_year, cutt.off.90days.scaled),
ep_stroke_t_90days = ifelse(ep_stroke_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_stroke_t_3years * avg_days_in_year, cutt.off.90days.scaled),
ep_coronary_t_90days = ifelse(ep_coronary_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_coronary_t_3years * avg_days_in_year, cutt.off.90days.scaled),
ep_cvdeath_t_90days = ifelse(ep_cvdeath_t_3years * avg_days_in_year <= cutt.off.90days.scaled,
ep_cvdeath_t_3years * avg_days_in_year, cutt.off.90days.scaled)
)
attach(AERNASE.clin.hdac9)
AERNASE.clin.hdac9[,"epmajor.30days"] <- AERNASE.clin.hdac9$epmajor.3years
AERNASE.clin.hdac9$epmajor.30days[epmajor.3years == 1 & ep_major_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epstroke.30days"] <- AERNASE.clin.hdac9$epstroke.3years
AERNASE.clin.hdac9$epstroke.30days[epstroke.3years == 1 & ep_stroke_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epcoronary.30days"] <- AERNASE.clin.hdac9$epcoronary.3years
AERNASE.clin.hdac9$epcoronary.30days[epcoronary.3years == 1 & ep_coronary_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epcvdeath.30days"] <- AERNASE.clin.hdac9$epcvdeath.3years
AERNASE.clin.hdac9$epcvdeath.30days[epcvdeath.3years == 1 & ep_cvdeath_t_3years > cutt.off.30days] <- 0
AERNASE.clin.hdac9[,"epmajor.90days"] <- AERNASE.clin.hdac9$epmajor.3years
AERNASE.clin.hdac9$epmajor.90days[epmajor.3years == 1 & ep_major_t_3years > cutt.off.90days] <- 0
AERNASE.clin.hdac9[,"epstroke.90days"] <- AERNASE.clin.hdac9$epstroke.3years
AERNASE.clin.hdac9$epstroke.90days[epstroke.3years == 1 & ep_stroke_t_3years > cutt.off.90days] <- 0
AERNASE.clin.hdac9[,"epcoronary.90days"] <- AERNASE.clin.hdac9$epcoronary.3years
AERNASE.clin.hdac9$epcoronary.90days[epcoronary.3years == 1 & ep_coronary_t_3years > cutt.off.90days] <- 0
AERNASE.clin.hdac9[,"epcvdeath.90days"] <- AERNASE.clin.hdac9$epcvdeath.3years
AERNASE.clin.hdac9$epcvdeath.90days[epcvdeath.3years == 1 & ep_cvdeath_t_3years > cutt.off.90days] <- 0
detach(AERNASE.clin.hdac9)
AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9, select = c("STUDY_NUMBER", "Age", "Gender", "Hospital", "Artery_summary",
"epmajor.3years", "epstroke.3years", "epcoronary.3years", "epcvdeath.3years",
"epmajor.30days", "epstroke.30days", "epcoronary.30days", "epcvdeath.30days",
"epmajor.90days", "epstroke.90days", "epcoronary.90days", "epcvdeath.90days"))
require(labelled)
AERNASE.clin.hdac9.temp$Gender <- to_factor(AERNASE.clin.hdac9.temp$Gender)
AERNASE.clin.hdac9.temp$Hospital <- to_factor(AERNASE.clin.hdac9.temp$Hospital)
AERNASE.clin.hdac9.temp$Artery_summary <- to_factor(AERNASE.clin.hdac9.temp$Artery_summary)
DT::datatable(AERNASE.clin.hdac9.temp[1:10,], caption = "Excerpt of the whole AERNASE.clin.hdac9.", rownames = FALSE)
rm(AERNASE.clin.hdac9.temp)First we do some sanity checks and inventory the time-to-event and event variables.
# Reference: https://bioconductor.org/packages/devel/bioc/vignettes/MultiAssayExperiment/inst/doc/QuickStartMultiAssay.html
# If you want to suppress warnings and messages when installing/loading packages
# suppressPackageStartupMessages({})
install.packages.auto("survival")
install.packages.auto("survminer")
install.packages.auto("Hmisc")
cat("* Creating function to summarize Cox regression and prepare container for results.")* Creating function to summarize Cox regression and prepare container for results.# Function to get summary statistics from Cox regression model
COX.STAT <- function(coxfit, DATASET, OUTCOME, target_of_interest){
cat("Summarizing Cox regression results for '", target_of_interest ,"' and its association to '",OUTCOME,"' in '",DATASET,"'.\n")
if (nrow(summary(coxfit)$coefficients) == 1) {
output = c(target_of_interest, rep(NA,8))
cat("Model not fitted; probably singular.\n")
}else {
cat("Collecting data.\n\n")
cox.sum <- summary(coxfit)
cox.effectsize = cox.sum$coefficients[1,1]
cox.SE = cox.sum$coefficients[1,3]
cox.HReffect = cox.sum$coefficients[1,2]
cox.CI_low = exp(cox.effectsize - 1.96 * cox.SE)
cox.CI_up = exp(cox.effectsize + 1.96 * cox.SE)
cox.zvalue = cox.sum$coefficients[1,4]
cox.pvalue = cox.sum$coefficients[1,5]
cox.sample_size = cox.sum$n
cox.nevents = cox.sum$nevent
output = c(DATASET, OUTCOME, target_of_interest, cox.effectsize, cox.SE, cox.HReffect, cox.CI_low, cox.CI_up, cox.zvalue, cox.pvalue, cox.sample_size, cox.nevents)
cat("We have collected the following:\n")
cat("Dataset used..............:", DATASET, "\n")
cat("Outcome analyzed..........:", OUTCOME, "\n")
cat("Protein...................:", target_of_interest, "\n")
cat("Effect size...............:", round(cox.effectsize, 6), "\n")
cat("Standard error............:", round(cox.SE, 6), "\n")
cat("Odds ratio (effect size)..:", round(cox.HReffect, 3), "\n")
cat("Lower 95% CI..............:", round(cox.CI_low, 3), "\n")
cat("Upper 95% CI..............:", round(cox.CI_up, 3), "\n")
cat("T-value...................:", round(cox.zvalue, 6), "\n")
cat("P-value...................:", signif(cox.pvalue, 8), "\n")
cat("Sample size in model......:", cox.sample_size, "\n")
cat("Number of events..........:", cox.nevents, "\n")
}
return(output)
print(output)
}
times = c("ep_major_t_3years",
"ep_stroke_t_3years", "ep_coronary_t_3years", "ep_cvdeath_t_3years")
endpoints = c("epmajor.3years",
"epstroke.3years", "epcoronary.3years", "epcvdeath.3years")
cat("* Check the cases per event type - for sanity.")* Check the cases per event type - for sanity.for (events in endpoints){
require(labelled)
print(paste0("Printing the summary of: ",events))
# print(summary(AERNASE.clin.hdac9[,events]))
print(table(AERNASE.clin.hdac9[,events]))
}[1] "Printing the summary of: epmajor.3years"
0 1
540 80
[1] "Printing the summary of: epstroke.3years"
0 1
573 47
[1] "Printing the summary of: epcoronary.3years"
0 1
572 48
[1] "Printing the summary of: epcvdeath.3years"
0 1
592 28 cat("* Check distribution of events over time - for sanity.")* Check distribution of events over time - for sanity.for (eventtimes in times){
print(paste0("Printing the summary of: ",eventtimes))
print(summary(AERNASE.clin.hdac9[,eventtimes]))
}[1] "Printing the summary of: ep_major_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.914 3.000 2.665 3.000 3.000 3
[1] "Printing the summary of: ep_stroke_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.929 3.000 2.708 3.000 3.000 3
[1] "Printing the summary of: ep_coronary_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.923 3.000 2.738 3.000 3.000 3
[1] "Printing the summary of: ep_cvdeath_t_3years"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0137 2.9500 3.0000 2.8246 3.0000 3.0000 3 for (eventtime in times){
print(paste0("Printing the distribution of: ",eventtime))
p <- gghistogram(AERNASE.clin.hdac9, x = eventtime, y = "..count..",
main = eventtime, bins = 15,
xlab = "year", color = uithof_color[16], fill = uithof_color[16], ggtheme = theme_minimal())
print(p)
ggsave(file = paste0(QC_loc, "/",Today,".AERNASE.clin.hdac9.EventDistributionPerYear.",eventtime,".pdf"), plot = last_plot())
}[1] "Printing the distribution of: ep_major_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_stroke_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_coronary_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_cvdeath_t_3years"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).
times30 = c("ep_major_t_30days",
"ep_stroke_t_30days", "ep_coronary_t_30days", "ep_cvdeath_t_30days")
endpoints30 = c("epmajor.30days",
"epstroke.30days", "epcoronary.30days", "epcvdeath.30days")
cat("* Check the cases per event type - for sanity.")* Check the cases per event type - for sanity.for (events in endpoints30){
print(paste0("Printing the summary of: ",events))
# print(summary(AERNASE.clin.hdac9[,events]))
print(table(AERNASE.clin.hdac9[,events]))
}[1] "Printing the summary of: epmajor.30days"
0 1
599 21
[1] "Printing the summary of: epstroke.30days"
0 1
602 18
[1] "Printing the summary of: epcoronary.30days"
0 1
614 6
[1] "Printing the summary of: epcvdeath.30days"
0 1
618 2 cat("* Check distribution of events over time - for sanity.")* Check distribution of events over time - for sanity.for (eventtimes in times30){
print(paste0("Printing the summary of: ",eventtimes))
print(summary(AERNASE.clin.hdac9[,eventtimes]))
}[1] "Printing the summary of: ep_major_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 30.00 30.00 29.05 30.00 30.00 3
[1] "Printing the summary of: ep_stroke_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 30.00 30.00 29.19 30.00 30.00 3
[1] "Printing the summary of: ep_coronary_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 30.00 30.00 29.67 30.00 30.00 3
[1] "Printing the summary of: ep_cvdeath_t_30days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
5.003 30.000 30.000 29.937 30.000 30.000 3 for (eventtime in times30){
print(paste0("Printing the distribution of: ",eventtime))
p <- gghistogram(AERNASE.clin.hdac9, x = eventtime, y = "..count..",
main = eventtime, bins = 15,
xlab = "days", color = uithof_color[16], fill = uithof_color[16], ggtheme = theme_minimal())
print(p)
ggsave(file = paste0(QC_loc, "/",Today,".AERNASE.clin.hdac9.EventDistributionPer30Days.",eventtime,".pdf"), plot = last_plot())
}[1] "Printing the distribution of: ep_major_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_stroke_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_coronary_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_cvdeath_t_30days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).
times90 = c("ep_major_t_90days",
"ep_stroke_t_90days", "ep_coronary_t_90days", "ep_cvdeath_t_90days")
endpoints90 = c("epmajor.90days",
"epstroke.90days", "epcoronary.90days", "epcvdeath.90days")
cat("* Check the cases per event type - for sanity.")* Check the cases per event type - for sanity.for (events in endpoints90){
print(paste0("Printing the summary of: ",events))
# print(summary(AERNASE.clin.hdac9[,events]))
print(table(AERNASE.clin.hdac9[,events]))
}[1] "Printing the summary of: epmajor.90days"
0 1
596 24
[1] "Printing the summary of: epstroke.90days"
0 1
600 20
[1] "Printing the summary of: epcoronary.90days"
0 1
612 8
[1] "Printing the summary of: epcvdeath.90days"
0 1
617 3 cat("* Check distribution of events over time - for sanity.")* Check distribution of events over time - for sanity.for (eventtimes in times90){
print(paste0("Printing the summary of: ",eventtimes))
print(summary(AERNASE.clin.hdac9[,eventtimes]))
}[1] "Printing the summary of: ep_major_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 90.00 90.00 86.79 90.00 90.00 3
[1] "Printing the summary of: ep_stroke_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 90.00 90.00 87.25 90.00 90.00 3
[1] "Printing the summary of: ep_coronary_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 90.00 90.00 88.75 90.00 90.00 3
[1] "Printing the summary of: ep_cvdeath_t_90days"
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
5.003 90.000 90.000 89.674 90.000 90.000 3 for (eventtime in times90){
print(paste0("Printing the distribution of: ",eventtime))
p <- gghistogram(AERNASE.clin.hdac9, x = eventtime, y = "..count..",
main = eventtime, bins = 15,
xlab = "days", color = uithof_color[16], fill = uithof_color[16], ggtheme = theme_minimal())
print(p)
ggsave(file = paste0(QC_loc, "/",Today,".AERNASE.clin.hdac9.EventDistributionPer90Days.",eventtime,".pdf"), plot = last_plot())
}[1] "Printing the distribution of: ep_major_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_stroke_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_coronary_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).[1] "Printing the distribution of: ep_cvdeath_t_90days"Warning: Removed 3 rows containing non-finite values (stat_bin).
Saving 18 x 12 in image
Warning: Removed 3 rows containing non-finite values (stat_bin).
Let’s perform the actual Cox-regressions. We will apply a couple of models:
# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times)){
eptime = times[i]
ep = endpoints[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [years]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 1 (Simple model)
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age+Gender + ORdate_year, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age+Gender + ORdate_year, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
# ylim = c(0.2, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [years]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL1.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}* Analyzing the effect of plaque target-of-interest on [epmajor.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12,13,14 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11,12,13,16 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 620, number of events= 80
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -3.351e-01 7.153e-01 2.322e-01 -1.443 0.1491
Age 3.283e-02 1.033e+00 1.377e-02 2.384 0.0171 *
Gendermale 7.816e-01 2.185e+00 3.270e-01 2.390 0.0168 *
ORdate_year2002 1.627e+01 1.165e+07 8.569e+03 0.002 0.9985
ORdate_year2003 1.513e+01 3.712e+06 8.569e+03 0.002 0.9986
ORdate_year2004 1.549e+01 5.320e+06 8.569e+03 0.002 0.9986
ORdate_year2005 1.582e+01 7.424e+06 8.569e+03 0.002 0.9985
ORdate_year2006 1.582e+01 7.416e+06 8.569e+03 0.002 0.9985
ORdate_year2007 1.497e+01 3.188e+06 8.569e+03 0.002 0.9986
ORdate_year2008 1.600e+01 8.864e+06 8.569e+03 0.002 0.9985
ORdate_year2009 1.533e+01 4.547e+06 8.569e+03 0.002 0.9986
ORdate_year2010 1.567e+01 6.404e+06 8.569e+03 0.002 0.9985
ORdate_year2011 1.475e+01 2.543e+06 8.569e+03 0.002 0.9986
ORdate_year2012 1.542e+01 4.972e+06 8.569e+03 0.002 0.9986
ORdate_year2013 -8.764e-01 4.163e-01 9.370e+03 0.000 0.9999
ORdate_year2014 -5.470e-01 5.787e-01 9.846e+03 0.000 1.0000
ORdate_year2015 1.640e-01 1.178e+00 9.788e+03 0.000 1.0000
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 7.153e-01 1.398e+00 0.4537 1.128
Age 1.033e+00 9.677e-01 1.0059 1.062
Gendermale 2.185e+00 4.577e-01 1.1511 4.148
ORdate_year2002 1.165e+07 8.586e-08 0.0000 Inf
ORdate_year2003 3.712e+06 2.694e-07 0.0000 Inf
ORdate_year2004 5.320e+06 1.880e-07 0.0000 Inf
ORdate_year2005 7.424e+06 1.347e-07 0.0000 Inf
ORdate_year2006 7.416e+06 1.348e-07 0.0000 Inf
ORdate_year2007 3.188e+06 3.136e-07 0.0000 Inf
ORdate_year2008 8.864e+06 1.128e-07 0.0000 Inf
ORdate_year2009 4.547e+06 2.199e-07 0.0000 Inf
ORdate_year2010 6.404e+06 1.561e-07 0.0000 Inf
ORdate_year2011 2.543e+06 3.933e-07 0.0000 Inf
ORdate_year2012 4.972e+06 2.011e-07 0.0000 Inf
ORdate_year2013 4.163e-01 2.402e+00 0.0000 Inf
ORdate_year2014 5.787e-01 1.728e+00 0.0000 Inf
ORdate_year2015 1.178e+00 8.488e-01 0.0000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Concordance= 0.656 (se = 0.029 )
Likelihood ratio test= 28.81 on 17 df, p=0.04
Wald test = 23.66 on 17 df, p=0.1
Score (logrank) test = 26.46 on 17 df, p=0.07
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epmajor.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epmajor.3years
Protein...................: HDAC9
Effect size...............: -0.335091
Standard error............: 0.232242
Odds ratio (effect size)..: 0.715
Lower 95% CI..............: 0.454
Upper 95% CI..............: 1.128
T-value...................: -1.442853
P-value...................: 0.1490618
Sample size in model......: 620
Number of events..........: 80
* Analyzing the effect of plaque target-of-interest on [epstroke.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12,13,14 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11,12,13 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 620, number of events= 47
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -3.870e-01 6.791e-01 3.038e-01 -1.274 0.203
Age 2.669e-02 1.027e+00 1.780e-02 1.499 0.134
Gendermale 3.709e-01 1.449e+00 3.741e-01 0.991 0.322
ORdate_year2002 1.657e+01 1.571e+07 1.150e+04 0.001 0.999
ORdate_year2003 1.557e+01 5.754e+06 1.150e+04 0.001 0.999
ORdate_year2004 1.575e+01 6.918e+06 1.150e+04 0.001 0.999
ORdate_year2005 1.625e+01 1.136e+07 1.150e+04 0.001 0.999
ORdate_year2006 1.650e+01 1.465e+07 1.150e+04 0.001 0.999
ORdate_year2007 1.501e+01 3.297e+06 1.150e+04 0.001 0.999
ORdate_year2008 1.637e+01 1.292e+07 1.150e+04 0.001 0.999
ORdate_year2009 1.506e+01 3.460e+06 1.150e+04 0.001 0.999
ORdate_year2010 1.641e+01 1.343e+07 1.150e+04 0.001 0.999
ORdate_year2011 1.502e+01 3.329e+06 1.150e+04 0.001 0.999
ORdate_year2012 1.602e+01 9.024e+06 1.150e+04 0.001 0.999
ORdate_year2013 -4.458e-01 6.403e-01 1.258e+04 0.000 1.000
ORdate_year2014 -1.032e-01 9.019e-01 1.321e+04 0.000 1.000
ORdate_year2015 3.256e-01 1.385e+00 1.324e+04 0.000 1.000
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 6.791e-01 1.473e+00 0.3744 1.232
Age 1.027e+00 9.737e-01 0.9918 1.064
Gendermale 1.449e+00 6.901e-01 0.6960 3.017
ORdate_year2002 1.571e+07 6.366e-08 0.0000 Inf
ORdate_year2003 5.754e+06 1.738e-07 0.0000 Inf
ORdate_year2004 6.918e+06 1.446e-07 0.0000 Inf
ORdate_year2005 1.136e+07 8.801e-08 0.0000 Inf
ORdate_year2006 1.465e+07 6.828e-08 0.0000 Inf
ORdate_year2007 3.297e+06 3.033e-07 0.0000 Inf
ORdate_year2008 1.292e+07 7.737e-08 0.0000 Inf
ORdate_year2009 3.460e+06 2.891e-07 0.0000 Inf
ORdate_year2010 1.343e+07 7.444e-08 0.0000 Inf
ORdate_year2011 3.329e+06 3.004e-07 0.0000 Inf
ORdate_year2012 9.024e+06 1.108e-07 0.0000 Inf
ORdate_year2013 6.403e-01 1.562e+00 0.0000 Inf
ORdate_year2014 9.019e-01 1.109e+00 0.0000 Inf
ORdate_year2015 1.385e+00 7.221e-01 0.0000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Concordance= 0.676 (se = 0.035 )
Likelihood ratio test= 19.37 on 17 df, p=0.3
Wald test = 13.33 on 17 df, p=0.7
Score (logrank) test = 17.81 on 17 df, p=0.4
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epstroke.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epstroke.3years
Protein...................: HDAC9
Effect size...............: -0.38697
Standard error............: 0.303845
Odds ratio (effect size)..: 0.679
Lower 95% CI..............: 0.374
Upper 95% CI..............: 1.232
T-value...................: -1.273579
P-value...................: 0.2028128
Sample size in model......: 620
Number of events..........: 47
* Analyzing the effect of plaque target-of-interest on [epcoronary.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12,13,14 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11,12,13 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 620, number of events= 48
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -3.056e-01 7.367e-01 2.966e-01 -1.030 0.3029
Age 1.465e-03 1.001e+00 1.721e-02 0.085 0.9321
Gendermale 9.151e-01 2.497e+00 4.391e-01 2.084 0.0372 *
ORdate_year2002 1.639e+01 1.313e+07 1.115e+04 0.001 0.9988
ORdate_year2003 1.578e+01 7.114e+06 1.115e+04 0.001 0.9989
ORdate_year2004 1.496e+01 3.140e+06 1.115e+04 0.001 0.9989
ORdate_year2005 1.571e+01 6.629e+06 1.115e+04 0.001 0.9989
ORdate_year2006 1.601e+01 8.936e+06 1.115e+04 0.001 0.9989
ORdate_year2007 1.553e+01 5.546e+06 1.115e+04 0.001 0.9989
ORdate_year2008 1.562e+01 6.058e+06 1.115e+04 0.001 0.9989
ORdate_year2009 1.533e+01 4.545e+06 1.115e+04 0.001 0.9989
ORdate_year2010 1.454e+01 2.063e+06 1.115e+04 0.001 0.9990
ORdate_year2011 1.464e+01 2.285e+06 1.115e+04 0.001 0.9990
ORdate_year2012 1.484e+01 2.782e+06 1.115e+04 0.001 0.9989
ORdate_year2013 -9.184e-01 3.991e-01 1.222e+04 0.000 0.9999
ORdate_year2014 -7.102e-01 4.915e-01 1.285e+04 0.000 1.0000
ORdate_year2015 -1.174e-01 8.892e-01 1.280e+04 0.000 1.0000
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 7.367e-01 1.357e+00 0.4119 1.318
Age 1.001e+00 9.985e-01 0.9682 1.036
Gendermale 2.497e+00 4.005e-01 1.0560 5.905
ORdate_year2002 1.313e+07 7.614e-08 0.0000 Inf
ORdate_year2003 7.114e+06 1.406e-07 0.0000 Inf
ORdate_year2004 3.140e+06 3.185e-07 0.0000 Inf
ORdate_year2005 6.629e+06 1.508e-07 0.0000 Inf
ORdate_year2006 8.936e+06 1.119e-07 0.0000 Inf
ORdate_year2007 5.546e+06 1.803e-07 0.0000 Inf
ORdate_year2008 6.058e+06 1.651e-07 0.0000 Inf
ORdate_year2009 4.545e+06 2.200e-07 0.0000 Inf
ORdate_year2010 2.063e+06 4.848e-07 0.0000 Inf
ORdate_year2011 2.285e+06 4.376e-07 0.0000 Inf
ORdate_year2012 2.782e+06 3.595e-07 0.0000 Inf
ORdate_year2013 3.991e-01 2.505e+00 0.0000 Inf
ORdate_year2014 4.915e-01 2.034e+00 0.0000 Inf
ORdate_year2015 8.892e-01 1.125e+00 0.0000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Concordance= 0.654 (se = 0.038 )
Likelihood ratio test= 17.21 on 17 df, p=0.4
Wald test = 13.48 on 17 df, p=0.7
Score (logrank) test = 15.71 on 17 df, p=0.5
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcoronary.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcoronary.3years
Protein...................: HDAC9
Effect size...............: -0.305628
Standard error............: 0.296646
Odds ratio (effect size)..: 0.737
Lower 95% CI..............: 0.412
Upper 95% CI..............: 1.318
T-value...................: -1.030276
P-value...................: 0.3028806
Sample size in model......: 620
Number of events..........: 48
* Analyzing the effect of plaque target-of-interest on [epcvdeath.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11 ; coefficient may be infinite. 
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year, data = TEMP.DF)
n= 620, number of events= 28
(3 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -7.408e-01 4.768e-01 4.175e-01 -1.774 0.07603 .
Age 7.826e-02 1.081e+00 2.610e-02 2.999 0.00271 **
Gendermale 6.706e-01 1.955e+00 5.434e-01 1.234 0.21717
ORdate_year2002 1.857e+01 1.158e+08 3.447e+04 0.001 0.99957
ORdate_year2003 1.695e+01 2.296e+07 3.447e+04 0.000 0.99961
ORdate_year2004 1.684e+01 2.054e+07 3.447e+04 0.000 0.99961
ORdate_year2005 1.756e+01 4.239e+07 3.447e+04 0.001 0.99959
ORdate_year2006 1.665e+01 1.701e+07 3.447e+04 0.000 0.99961
ORdate_year2007 1.720e+01 2.962e+07 3.447e+04 0.000 0.99960
ORdate_year2008 1.730e+01 3.257e+07 3.447e+04 0.001 0.99960
ORdate_year2009 1.758e+01 4.302e+07 3.447e+04 0.001 0.99959
ORdate_year2010 1.749e+01 3.927e+07 3.447e+04 0.001 0.99960
ORdate_year2011 -5.884e-01 5.552e-01 3.487e+04 0.000 0.99999
ORdate_year2012 -6.736e-01 5.099e-01 3.507e+04 0.000 0.99998
ORdate_year2013 -9.400e-01 3.906e-01 3.728e+04 0.000 0.99998
ORdate_year2014 -1.580e-01 8.539e-01 3.894e+04 0.000 1.00000
ORdate_year2015 7.661e-01 2.151e+00 3.907e+04 0.000 0.99998
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 4.768e-01 2.098e+00 0.2103 1.081
Age 1.081e+00 9.247e-01 1.0275 1.138
Gendermale 1.955e+00 5.114e-01 0.6740 5.673
ORdate_year2002 1.158e+08 8.639e-09 0.0000 Inf
ORdate_year2003 2.296e+07 4.356e-08 0.0000 Inf
ORdate_year2004 2.054e+07 4.868e-08 0.0000 Inf
ORdate_year2005 4.239e+07 2.359e-08 0.0000 Inf
ORdate_year2006 1.701e+07 5.880e-08 0.0000 Inf
ORdate_year2007 2.962e+07 3.376e-08 0.0000 Inf
ORdate_year2008 3.257e+07 3.070e-08 0.0000 Inf
ORdate_year2009 4.302e+07 2.325e-08 0.0000 Inf
ORdate_year2010 3.927e+07 2.546e-08 0.0000 Inf
ORdate_year2011 5.552e-01 1.801e+00 0.0000 Inf
ORdate_year2012 5.099e-01 1.961e+00 0.0000 Inf
ORdate_year2013 3.906e-01 2.560e+00 0.0000 Inf
ORdate_year2014 8.539e-01 1.171e+00 0.0000 Inf
ORdate_year2015 2.151e+00 4.648e-01 0.0000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Concordance= 0.792 (se = 0.037 )
Likelihood ratio test= 30.55 on 17 df, p=0.02
Wald test = 14.69 on 17 df, p=0.6
Score (logrank) test = 30.08 on 17 df, p=0.03
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcvdeath.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcvdeath.3years
Protein...................: HDAC9
Effect size...............: -0.740749
Standard error............: 0.417516
Odds ratio (effect size)..: 0.477
Lower 95% CI..............: 0.21
Upper 95% CI..............: 1.081
T-value...................: -1.77418
P-value...................: 0.07603341
Sample size in model......: 620
Number of events..........: 28 
cat("- Edit the column names...\n")- Edit the column names...colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")- Correct the variable types...COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")- Writing results to Excel-file...head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL1.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")- Removing intermediate files...rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times)){
eptime = times[i]
ep = endpoints[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
#as.integer(TEMP.DF[,ep] == "Excluded")
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [years]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 2 adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
# ylim = c(0.2, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [years]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL2.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}* Analyzing the effect of plaque target-of-interest on [epmajor.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12,13,14,35,36,37 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11,12,13,34,35,36 ; coefficient may be infinite.
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 535, number of events= 65
(88 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -1.070e-01 8.985e-01 2.602e-01 -0.411 0.6808
Age 4.732e-02 1.048e+00 1.856e-02 2.550 0.0108 *
Gendermale 8.189e-01 2.268e+00 3.879e-01 2.111 0.0348 *
ORdate_year2002 1.536e+01 4.678e+06 1.019e+04 0.002 0.9988
ORdate_year2003 1.429e+01 1.601e+06 1.019e+04 0.001 0.9989
ORdate_year2004 1.475e+01 2.550e+06 1.019e+04 0.001 0.9988
ORdate_year2005 1.523e+01 4.116e+06 1.019e+04 0.001 0.9988
ORdate_year2006 1.520e+01 4.007e+06 1.019e+04 0.001 0.9988
ORdate_year2007 1.422e+01 1.495e+06 1.019e+04 0.001 0.9989
ORdate_year2008 1.528e+01 4.311e+06 1.019e+04 0.001 0.9988
ORdate_year2009 1.418e+01 1.446e+06 1.019e+04 0.001 0.9989
ORdate_year2010 1.508e+01 3.550e+06 1.019e+04 0.001 0.9988
ORdate_year2011 1.440e+01 1.795e+06 1.019e+04 0.001 0.9989
ORdate_year2012 1.498e+01 3.205e+06 1.019e+04 0.001 0.9988
ORdate_year2013 -1.394e+00 2.480e-01 1.121e+04 0.000 0.9999
ORdate_year2014 -9.861e-01 3.730e-01 1.441e+04 0.000 0.9999
ORdate_year2015 -7.423e-01 4.760e-01 1.128e+04 0.000 0.9999
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
Hypertension.compositeno -6.272e-01 5.341e-01 5.432e-01 -1.155 0.2482
Hypertension.compositeyes NA NA 0.000e+00 NA NA
DiabetesStatusDiabetes 6.858e-01 1.985e+00 2.800e-01 2.449 0.0143 *
SmokerStatusEx-smoker -5.207e-01 5.941e-01 2.781e-01 -1.873 0.0611 .
SmokerStatusNever smoked -8.051e-01 4.470e-01 4.407e-01 -1.827 0.0677 .
Med.Statin.LLDno 2.396e-01 1.271e+00 2.979e-01 0.804 0.4213
Med.Statin.LLDyes NA NA 0.000e+00 NA NA
Med.all.antiplateletno 1.650e-01 1.179e+00 4.249e-01 0.388 0.6978
Med.all.antiplateletyes NA NA 0.000e+00 NA NA
GFR_MDRD -7.953e-03 9.921e-01 6.443e-03 -1.234 0.2170
BMI 5.476e-02 1.056e+00 3.163e-02 1.731 0.0834 .
MedHx_CVDNo -5.580e-01 5.724e-01 2.995e-01 -1.863 0.0624 .
stenose0-49% 4.153e-01 1.515e+00 1.156e+04 0.000 1.0000
stenose50-70% 1.610e+01 9.869e+06 1.019e+04 0.002 0.9987
stenose70-90% 1.663e+01 1.669e+07 1.019e+04 0.002 0.9987
stenose90-99% 1.671e+01 1.802e+07 1.019e+04 0.002 0.9987
stenose100% (Occlusion) 1.043e+00 2.837e+00 1.113e+04 0.000 0.9999
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% NA NA 0.000e+00 NA NA
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 8.985e-01 1.113e+00 0.5396 1.496
Age 1.048e+00 9.538e-01 1.0110 1.087
Gendermale 2.268e+00 4.409e-01 1.0603 4.851
ORdate_year2002 4.678e+06 2.138e-07 0.0000 Inf
ORdate_year2003 1.601e+06 6.247e-07 0.0000 Inf
ORdate_year2004 2.550e+06 3.922e-07 0.0000 Inf
ORdate_year2005 4.116e+06 2.430e-07 0.0000 Inf
ORdate_year2006 4.007e+06 2.495e-07 0.0000 Inf
ORdate_year2007 1.495e+06 6.687e-07 0.0000 Inf
ORdate_year2008 4.311e+06 2.320e-07 0.0000 Inf
ORdate_year2009 1.446e+06 6.915e-07 0.0000 Inf
ORdate_year2010 3.550e+06 2.817e-07 0.0000 Inf
ORdate_year2011 1.795e+06 5.571e-07 0.0000 Inf
ORdate_year2012 3.205e+06 3.120e-07 0.0000 Inf
ORdate_year2013 2.480e-01 4.033e+00 0.0000 Inf
ORdate_year2014 3.730e-01 2.681e+00 0.0000 Inf
ORdate_year2015 4.760e-01 2.101e+00 0.0000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Hypertension.compositeno 5.341e-01 1.872e+00 0.1842 1.549
Hypertension.compositeyes NA NA NA NA
DiabetesStatusDiabetes 1.985e+00 5.037e-01 1.1468 3.437
SmokerStatusEx-smoker 5.941e-01 1.683e+00 0.3445 1.025
SmokerStatusNever smoked 4.470e-01 2.237e+00 0.1885 1.060
Med.Statin.LLDno 1.271e+00 7.870e-01 0.7087 2.279
Med.Statin.LLDyes NA NA NA NA
Med.all.antiplateletno 1.179e+00 8.479e-01 0.5128 2.712
Med.all.antiplateletyes NA NA NA NA
GFR_MDRD 9.921e-01 1.008e+00 0.9796 1.005
BMI 1.056e+00 9.467e-01 0.9928 1.124
MedHx_CVDNo 5.724e-01 1.747e+00 0.3183 1.029
stenose0-49% 1.515e+00 6.601e-01 0.0000 Inf
stenose50-70% 9.869e+06 1.013e-07 0.0000 Inf
stenose70-90% 1.669e+07 5.991e-08 0.0000 Inf
stenose90-99% 1.802e+07 5.550e-08 0.0000 Inf
stenose100% (Occlusion) 2.837e+00 3.525e-01 0.0000 Inf
stenoseNA NA NA NA NA
stenose50-99% NA NA NA NA
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.764 (se = 0.025 )
Likelihood ratio test= 55.11 on 31 df, p=0.005
Wald test = 38.75 on 31 df, p=0.2
Score (logrank) test = 53.61 on 31 df, p=0.007
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epmajor.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epmajor.3years
Protein...................: HDAC9
Effect size...............: -0.10702
Standard error............: 0.260153
Odds ratio (effect size)..: 0.899
Lower 95% CI..............: 0.54
Upper 95% CI..............: 1.496
T-value...................: -0.411374
P-value...................: 0.6807986
Sample size in model......: 535
Number of events..........: 65
* Analyzing the effect of plaque target-of-interest on [epstroke.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,12,13,14,35,36,37 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,11,12,13,34,35,36 ; coefficient may be infinite.
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 535, number of events= 34
(88 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -1.405e-01 8.689e-01 3.612e-01 -0.389 0.6973
Age 5.009e-02 1.051e+00 2.480e-02 2.020 0.0434 *
Gendermale 4.427e-01 1.557e+00 4.682e-01 0.945 0.3444
ORdate_year2002 1.745e+01 3.778e+07 4.162e+04 0.000 0.9997
ORdate_year2003 1.639e+01 1.318e+07 4.162e+04 0.000 0.9997
ORdate_year2004 1.676e+01 1.897e+07 4.162e+04 0.000 0.9997
ORdate_year2005 1.730e+01 3.260e+07 4.162e+04 0.000 0.9997
ORdate_year2006 1.782e+01 5.483e+07 4.162e+04 0.000 0.9997
ORdate_year2007 1.574e+01 6.883e+06 4.162e+04 0.000 0.9997
ORdate_year2008 1.732e+01 3.332e+07 4.162e+04 0.000 0.9997
ORdate_year2009 -1.313e+00 2.691e-01 4.190e+04 0.000 1.0000
ORdate_year2010 1.758e+01 4.300e+07 4.162e+04 0.000 0.9997
ORdate_year2011 1.647e+01 1.420e+07 4.162e+04 0.000 0.9997
ORdate_year2012 1.725e+01 3.089e+07 4.162e+04 0.000 0.9997
ORdate_year2013 -1.102e+00 3.321e-01 4.614e+04 0.000 1.0000
ORdate_year2014 -9.220e-01 3.977e-01 5.885e+04 0.000 1.0000
ORdate_year2015 -9.160e-01 4.001e-01 4.734e+04 0.000 1.0000
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
Hypertension.compositeno -1.841e-01 8.319e-01 6.509e-01 -0.283 0.7773
Hypertension.compositeyes NA NA 0.000e+00 NA NA
DiabetesStatusDiabetes 3.340e-01 1.397e+00 4.180e-01 0.799 0.4243
SmokerStatusEx-smoker -5.879e-01 5.555e-01 3.743e-01 -1.571 0.1162
SmokerStatusNever smoked -1.568e+00 2.085e-01 7.651e-01 -2.049 0.0404 *
Med.Statin.LLDno 2.917e-01 1.339e+00 4.067e-01 0.717 0.4733
Med.Statin.LLDyes NA NA 0.000e+00 NA NA
Med.all.antiplateletno -6.841e-02 9.339e-01 6.452e-01 -0.106 0.9156
Med.all.antiplateletyes NA NA 0.000e+00 NA NA
GFR_MDRD 6.356e-03 1.006e+00 8.677e-03 0.732 0.4639
BMI 9.627e-02 1.101e+00 4.095e-02 2.351 0.0187 *
MedHx_CVDNo -4.748e-01 6.220e-01 4.047e-01 -1.173 0.2406
stenose0-49% 1.588e+00 4.896e+00 4.579e+04 0.000 1.0000
stenose50-70% 1.952e+01 2.998e+08 4.161e+04 0.000 0.9996
stenose70-90% 1.942e+01 2.725e+08 4.161e+04 0.000 0.9996
stenose90-99% 1.951e+01 2.959e+08 4.161e+04 0.000 0.9996
stenose100% (Occlusion) 1.617e+00 5.039e+00 4.436e+04 0.000 1.0000
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% NA NA 0.000e+00 NA NA
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 8.689e-01 1.151e+00 0.42810 1.7636
Age 1.051e+00 9.511e-01 1.00148 1.1037
Gendermale 1.557e+00 6.423e-01 0.62185 3.8978
ORdate_year2002 3.778e+07 2.647e-08 0.00000 Inf
ORdate_year2003 1.318e+07 7.589e-08 0.00000 Inf
ORdate_year2004 1.897e+07 5.272e-08 0.00000 Inf
ORdate_year2005 3.260e+07 3.068e-08 0.00000 Inf
ORdate_year2006 5.483e+07 1.824e-08 0.00000 Inf
ORdate_year2007 6.883e+06 1.453e-07 0.00000 Inf
ORdate_year2008 3.332e+07 3.001e-08 0.00000 Inf
ORdate_year2009 2.691e-01 3.716e+00 0.00000 Inf
ORdate_year2010 4.300e+07 2.325e-08 0.00000 Inf
ORdate_year2011 1.420e+07 7.044e-08 0.00000 Inf
ORdate_year2012 3.089e+07 3.237e-08 0.00000 Inf
ORdate_year2013 3.321e-01 3.011e+00 0.00000 Inf
ORdate_year2014 3.977e-01 2.514e+00 0.00000 Inf
ORdate_year2015 4.001e-01 2.499e+00 0.00000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Hypertension.compositeno 8.319e-01 1.202e+00 0.23227 2.9793
Hypertension.compositeyes NA NA NA NA
DiabetesStatusDiabetes 1.397e+00 7.161e-01 0.61551 3.1686
SmokerStatusEx-smoker 5.555e-01 1.800e+00 0.26673 1.1568
SmokerStatusNever smoked 2.085e-01 4.797e+00 0.04654 0.9339
Med.Statin.LLDno 1.339e+00 7.470e-01 0.60326 2.9706
Med.Statin.LLDyes NA NA NA NA
Med.all.antiplateletno 9.339e-01 1.071e+00 0.26371 3.3071
Med.all.antiplateletyes NA NA NA NA
GFR_MDRD 1.006e+00 9.937e-01 0.98941 1.0236
BMI 1.101e+00 9.082e-01 1.01613 1.1931
MedHx_CVDNo 6.220e-01 1.608e+00 0.28141 1.3748
stenose0-49% 4.896e+00 2.043e-01 0.00000 Inf
stenose50-70% 2.998e+08 3.335e-09 0.00000 Inf
stenose70-90% 2.725e+08 3.669e-09 0.00000 Inf
stenose90-99% 2.959e+08 3.379e-09 0.00000 Inf
stenose100% (Occlusion) 5.039e+00 1.985e-01 0.00000 Inf
stenoseNA NA NA NA NA
stenose50-99% NA NA NA NA
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.795 (se = 0.031 )
Likelihood ratio test= 39.43 on 31 df, p=0.1
Wald test = 16.75 on 31 df, p=1
Score (logrank) test = 35.86 on 31 df, p=0.3
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epstroke.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epstroke.3years
Protein...................: HDAC9
Effect size...............: -0.140508
Standard error............: 0.361176
Odds ratio (effect size)..: 0.869
Lower 95% CI..............: 0.428
Upper 95% CI..............: 1.764
T-value...................: -0.38903
P-value...................: 0.6972542
Sample size in model......: 535
Number of events..........: 34
* Analyzing the effect of plaque target-of-interest on [epcoronary.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12,13,14,16,17,35,36,37,38 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11,12,13,35,36,37 ; coefficient may be infinite.
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 535, number of events= 43
(88 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -1.546e-01 8.567e-01 3.205e-01 -0.482 0.62953
Age -1.299e-03 9.987e-01 2.280e-02 -0.057 0.95459
Gendermale 1.273e+00 3.572e+00 5.389e-01 2.362 0.01816 *
ORdate_year2002 1.746e+01 3.822e+07 3.427e+04 0.001 0.99959
ORdate_year2003 1.707e+01 2.601e+07 3.427e+04 0.000 0.99960
ORdate_year2004 1.624e+01 1.132e+07 3.427e+04 0.000 0.99962
ORdate_year2005 1.722e+01 3.015e+07 3.427e+04 0.001 0.99960
ORdate_year2006 1.702e+01 2.464e+07 3.427e+04 0.000 0.99960
ORdate_year2007 1.684e+01 2.054e+07 3.427e+04 0.000 0.99961
ORdate_year2008 1.697e+01 2.344e+07 3.427e+04 0.000 0.99960
ORdate_year2009 1.670e+01 1.784e+07 3.427e+04 0.000 0.99961
ORdate_year2010 1.563e+01 6.123e+06 3.427e+04 0.000 0.99964
ORdate_year2011 1.602e+01 9.103e+06 3.427e+04 0.000 0.99963
ORdate_year2012 1.654e+01 1.525e+07 3.427e+04 0.000 0.99961
ORdate_year2013 -1.790e+00 1.669e-01 3.758e+04 0.000 0.99996
ORdate_year2014 -5.827e-01 5.584e-01 4.847e+04 0.000 0.99999
ORdate_year2015 -1.527e-01 8.584e-01 3.747e+04 0.000 1.00000
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
Hypertension.compositeno -1.679e+00 1.865e-01 1.029e+00 -1.631 0.10283
Hypertension.compositeyes NA NA 0.000e+00 NA NA
DiabetesStatusDiabetes 3.192e-01 1.376e+00 3.502e-01 0.912 0.36194
SmokerStatusEx-smoker -3.283e-01 7.201e-01 3.536e-01 -0.929 0.35308
SmokerStatusNever smoked 9.937e-02 1.104e+00 4.722e-01 0.210 0.83332
Med.Statin.LLDno -4.074e-01 6.654e-01 4.369e-01 -0.932 0.35113
Med.Statin.LLDyes NA NA 0.000e+00 NA NA
Med.all.antiplateletno 1.760e-02 1.018e+00 5.674e-01 0.031 0.97526
Med.all.antiplateletyes NA NA 0.000e+00 NA NA
GFR_MDRD -1.438e-02 9.857e-01 8.199e-03 -1.754 0.07940 .
BMI 1.769e-02 1.018e+00 4.199e-02 0.421 0.67352
MedHx_CVDNo -1.199e+00 3.016e-01 4.254e-01 -2.818 0.00484 **
stenose0-49% -1.731e+00 1.771e-01 3.919e+04 0.000 0.99996
stenose50-70% -6.636e-01 5.150e-01 3.469e+04 0.000 0.99998
stenose70-90% 1.726e+01 3.129e+07 3.427e+04 0.001 0.99960
stenose90-99% 1.708e+01 2.616e+07 3.427e+04 0.000 0.99960
stenose100% (Occlusion) -5.124e-01 5.990e-01 3.719e+04 0.000 0.99999
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% NA NA 0.000e+00 NA NA
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 8.567e-01 1.167e+00 0.4571 1.6057
Age 9.987e-01 1.001e+00 0.9551 1.0443
Gendermale 3.572e+00 2.800e-01 1.2421 10.2700
ORdate_year2002 3.822e+07 2.617e-08 0.0000 Inf
ORdate_year2003 2.601e+07 3.845e-08 0.0000 Inf
ORdate_year2004 1.132e+07 8.835e-08 0.0000 Inf
ORdate_year2005 3.015e+07 3.317e-08 0.0000 Inf
ORdate_year2006 2.464e+07 4.059e-08 0.0000 Inf
ORdate_year2007 2.054e+07 4.869e-08 0.0000 Inf
ORdate_year2008 2.344e+07 4.266e-08 0.0000 Inf
ORdate_year2009 1.784e+07 5.606e-08 0.0000 Inf
ORdate_year2010 6.123e+06 1.633e-07 0.0000 Inf
ORdate_year2011 9.103e+06 1.099e-07 0.0000 Inf
ORdate_year2012 1.525e+07 6.557e-08 0.0000 Inf
ORdate_year2013 1.669e-01 5.992e+00 0.0000 Inf
ORdate_year2014 5.584e-01 1.791e+00 0.0000 Inf
ORdate_year2015 8.584e-01 1.165e+00 0.0000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Hypertension.compositeno 1.865e-01 5.362e+00 0.0248 1.4026
Hypertension.compositeyes NA NA NA NA
DiabetesStatusDiabetes 1.376e+00 7.267e-01 0.6928 2.7333
SmokerStatusEx-smoker 7.201e-01 1.389e+00 0.3601 1.4400
SmokerStatusNever smoked 1.104e+00 9.054e-01 0.4378 2.7866
Med.Statin.LLDno 6.654e-01 1.503e+00 0.2826 1.5667
Med.Statin.LLDyes NA NA NA NA
Med.all.antiplateletno 1.018e+00 9.826e-01 0.3347 3.0945
Med.all.antiplateletyes NA NA NA NA
GFR_MDRD 9.857e-01 1.014e+00 0.9700 1.0017
BMI 1.018e+00 9.825e-01 0.9374 1.1051
MedHx_CVDNo 3.016e-01 3.316e+00 0.1310 0.6943
stenose0-49% 1.771e-01 5.646e+00 0.0000 Inf
stenose50-70% 5.150e-01 1.942e+00 0.0000 Inf
stenose70-90% 3.129e+07 3.196e-08 0.0000 Inf
stenose90-99% 2.616e+07 3.822e-08 0.0000 Inf
stenose100% (Occlusion) 5.990e-01 1.669e+00 0.0000 Inf
stenoseNA NA NA NA NA
stenose50-99% NA NA NA NA
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.778 (se = 0.027 )
Likelihood ratio test= 45.51 on 31 df, p=0.04
Wald test = 23.94 on 31 df, p=0.8
Score (logrank) test = 37.26 on 31 df, p=0.2
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcoronary.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcoronary.3years
Protein...................: HDAC9
Effect size...............: -0.15461
Standard error............: 0.320506
Odds ratio (effect size)..: 0.857
Lower 95% CI..............: 0.457
Upper 95% CI..............: 1.606
T-value...................: -0.482393
P-value...................: 0.6295269
Sample size in model......: 535
Number of events..........: 43
* Analyzing the effect of plaque target-of-interest on [epcvdeath.3years].
- creating temporary SE for this work.
- making a 'Surv' object and adding this to temporary dataframe.
- making strata of each of the plaque target-of-interest and start survival analysis.
> processing [HDAC9]; 1 out of 1 target-of-interest.
> cross tabulation of HDAC9-stratum.
[ 0, 22) [22,449]
319 304
> fitting the model for HDAC9-stratum.
> make a Kaplan-Meier-shizzle...
> perform the Cox-regression fashizzle and plot it...Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 4,5,6,7,8,9,10,11,12,13,14,15,16,17,22,34,35,36,37,38 ; coefficient may be infinite.
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 3,4,5,6,7,8,9,10,11,13,16,21,33,34,35,36,37 ; coefficient may be infinite.
Call:
coxph(formula = Surv(TEMP.DF[, eptime], event) ~ TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]] +
Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus +
SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD +
BMI + MedHx_CVD + stenose, data = TEMP.DF)
n= 535, number of events= 23
(88 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] -2.771e-01 7.580e-01 4.600e-01 -0.602 0.5469
Age 7.215e-02 1.075e+00 3.477e-02 2.075 0.0380 *
Gendermale 6.159e-01 1.851e+00 6.573e-01 0.937 0.3487
ORdate_year2002 1.820e+01 8.034e+07 7.481e+04 0.000 0.9998
ORdate_year2003 1.684e+01 2.065e+07 7.481e+04 0.000 0.9998
ORdate_year2004 1.699e+01 2.393e+07 7.481e+04 0.000 0.9998
ORdate_year2005 1.775e+01 5.101e+07 7.481e+04 0.000 0.9998
ORdate_year2006 1.673e+01 1.847e+07 7.481e+04 0.000 0.9998
ORdate_year2007 1.709e+01 2.642e+07 7.481e+04 0.000 0.9998
ORdate_year2008 1.726e+01 3.131e+07 7.481e+04 0.000 0.9998
ORdate_year2009 1.691e+01 2.214e+07 7.481e+04 0.000 0.9998
ORdate_year2010 1.766e+01 4.657e+07 7.481e+04 0.000 0.9998
ORdate_year2011 -1.550e+00 2.123e-01 7.545e+04 0.000 1.0000
ORdate_year2012 -1.213e+00 2.974e-01 7.591e+04 0.000 1.0000
ORdate_year2013 -1.558e+00 2.105e-01 8.214e+04 0.000 1.0000
ORdate_year2014 -1.108e+00 3.303e-01 1.058e+05 0.000 1.0000
ORdate_year2015 -1.084e+00 3.383e-01 7.997e+04 0.000 1.0000
ORdate_year2016 NA NA 0.000e+00 NA NA
ORdate_year2017 NA NA 0.000e+00 NA NA
ORdate_year2018 NA NA 0.000e+00 NA NA
ORdate_year2019 NA NA 0.000e+00 NA NA
Hypertension.compositeno -1.820e+01 1.248e-08 4.924e+03 -0.004 0.9971
Hypertension.compositeyes NA NA 0.000e+00 NA NA
DiabetesStatusDiabetes 9.437e-01 2.570e+00 4.835e-01 1.952 0.0510 .
SmokerStatusEx-smoker -6.808e-01 5.062e-01 4.963e-01 -1.372 0.1701
SmokerStatusNever smoked -1.101e+00 3.325e-01 8.341e-01 -1.320 0.1868
Med.Statin.LLDno 3.260e-01 1.385e+00 4.851e-01 0.672 0.5015
Med.Statin.LLDyes NA NA 0.000e+00 NA NA
Med.all.antiplateletno 1.089e+00 2.970e+00 5.963e-01 1.826 0.0679 .
Med.all.antiplateletyes NA NA 0.000e+00 NA NA
GFR_MDRD -2.453e-02 9.758e-01 1.207e-02 -2.032 0.0422 *
BMI 3.188e-02 1.032e+00 6.026e-02 0.529 0.5968
MedHx_CVDNo -5.370e-01 5.845e-01 5.478e-01 -0.980 0.3269
stenose0-49% -3.568e+01 3.197e-16 8.202e+04 0.000 0.9997
stenose50-70% -1.702e+01 4.059e-08 7.607e+04 0.000 0.9998
stenose70-90% -1.733e+01 2.984e-08 7.607e+04 0.000 0.9998
stenose90-99% -1.703e+01 4.016e-08 7.607e+04 0.000 0.9998
stenose100% (Occlusion) -3.299e+01 4.695e-15 7.795e+04 0.000 0.9997
stenoseNA NA NA 0.000e+00 NA NA
stenose50-99% NA NA 0.000e+00 NA NA
stenose70-99% NA NA 0.000e+00 NA NA
stenose99 NA NA 0.000e+00 NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
TEMP.DF[[TRAITS.TARGET.RANK[target_of_interest]]][22,449] 7.580e-01 1.319e+00 0.30766 1.8673
Age 1.075e+00 9.304e-01 1.00401 1.1506
Gendermale 1.851e+00 5.401e-01 0.51053 6.7140
ORdate_year2002 8.034e+07 1.245e-08 0.00000 Inf
ORdate_year2003 2.065e+07 4.844e-08 0.00000 Inf
ORdate_year2004 2.393e+07 4.179e-08 0.00000 Inf
ORdate_year2005 5.101e+07 1.960e-08 0.00000 Inf
ORdate_year2006 1.847e+07 5.416e-08 0.00000 Inf
ORdate_year2007 2.642e+07 3.785e-08 0.00000 Inf
ORdate_year2008 3.131e+07 3.194e-08 0.00000 Inf
ORdate_year2009 2.214e+07 4.516e-08 0.00000 Inf
ORdate_year2010 4.657e+07 2.147e-08 0.00000 Inf
ORdate_year2011 2.123e-01 4.710e+00 0.00000 Inf
ORdate_year2012 2.974e-01 3.363e+00 0.00000 Inf
ORdate_year2013 2.105e-01 4.750e+00 0.00000 Inf
ORdate_year2014 3.303e-01 3.027e+00 0.00000 Inf
ORdate_year2015 3.383e-01 2.956e+00 0.00000 Inf
ORdate_year2016 NA NA NA NA
ORdate_year2017 NA NA NA NA
ORdate_year2018 NA NA NA NA
ORdate_year2019 NA NA NA NA
Hypertension.compositeno 1.248e-08 8.016e+07 0.00000 Inf
Hypertension.compositeyes NA NA NA NA
DiabetesStatusDiabetes 2.570e+00 3.892e-01 0.99604 6.6287
SmokerStatusEx-smoker 5.062e-01 1.975e+00 0.19137 1.3390
SmokerStatusNever smoked 3.325e-01 3.007e+00 0.06484 1.7051
Med.Statin.LLDno 1.385e+00 7.218e-01 0.53540 3.5853
Med.Statin.LLDyes NA NA NA NA
Med.all.antiplateletno 2.970e+00 3.367e-01 0.92302 9.5590
Med.all.antiplateletyes NA NA NA NA
GFR_MDRD 9.758e-01 1.025e+00 0.95295 0.9991
BMI 1.032e+00 9.686e-01 0.91739 1.1618
MedHx_CVDNo 5.845e-01 1.711e+00 0.19978 1.7102
stenose0-49% 3.197e-16 3.128e+15 0.00000 Inf
stenose50-70% 4.059e-08 2.464e+07 0.00000 Inf
stenose70-90% 2.984e-08 3.351e+07 0.00000 Inf
stenose90-99% 4.016e-08 2.490e+07 0.00000 Inf
stenose100% (Occlusion) 4.695e-15 2.130e+14 0.00000 Inf
stenoseNA NA NA NA NA
stenose50-99% NA NA NA NA
stenose70-99% NA NA NA NA
stenose99 NA NA NA NA
Concordance= 0.852 (se = 0.033 )
Likelihood ratio test= 42.92 on 31 df, p=0.08
Wald test = 21.37 on 31 df, p=0.9
Score (logrank) test = 39.69 on 31 df, p=0.1
> writing the Cox-regression fashizzle to Excel...
Summarizing Cox regression results for ' HDAC9 ' and its association to ' epcvdeath.3years ' in ' AERNASE.clin.hdac9 '.
Collecting data.
We have collected the following:
Dataset used..............: AERNASE.clin.hdac9
Outcome analyzed..........: epcvdeath.3years
Protein...................: HDAC9
Effect size...............: -0.277124
Standard error............: 0.46002
Odds ratio (effect size)..: 0.758
Lower 95% CI..............: 0.308
Upper 95% CI..............: 1.867
T-value...................: -0.602416
P-value...................: 0.546897
Sample size in model......: 535
Number of events..........: 23 
cat("- Edit the column names...\n")- Edit the column names...colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")- Correct the variable types...COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")- Writing results to Excel-file...head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL2.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")- Removing intermediate files...rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times30)){
eptime = times30[i]
ep = endpoints30[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".30days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 1 (Simple model)
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age+Gender + ORdate_year, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age+Gender + ORdate_year, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL1.30days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
library(openxlsx)
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL1.30days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colnames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times30)){
eptime = times30[i]
ep = endpoints30[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
#as.integer(TEMP.DF[,ep] == "Excluded")
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".30days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 2 adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL2.30days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL2.30days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times90)){
eptime = times90[i]
ep = endpoints90[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".90days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 1 (Simple model)
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age+Gender + ORdate_year, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age+Gender + ORdate_year, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#595A5C", "#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL1.90days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
library(openxlsx)
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL1.90days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)# Set up a dataframe to receive results
COX.results <- data.frame(matrix(NA, ncol = 12, nrow = 0))
# Looping over each target_of_interest/endpoint/time combination
for (i in 1:length(times90)){
eptime = times90[i]
ep = endpoints90[i]
cat(paste0("* Analyzing the effect of plaque target-of-interest on [",ep,"].\n"))
cat(" - creating temporary SE for this work.\n")
TEMP.DF = as.data.frame(AERNASE.clin.hdac9)
cat(" - making a 'Surv' object and adding this to temporary dataframe.\n")
TEMP.DF$event <- as.integer(TEMP.DF[,ep])
#as.integer(TEMP.DF[,ep] == "Excluded")
TEMP.DF$y <- Surv(time = TEMP.DF[,eptime], event = TEMP.DF$event)
cat(" - making strata of each of the plaque target-of-interest and start survival analysis.\n")
for (target_of_interest in 1:length(TRAITS.TARGET.RANK)){
cat(paste0(" > processing [",TRAITS.TARGET.RANK[target_of_interest],"]; ",target_of_interest," out of ",length(TRAITS.TARGET.RANK)," target-of-interest.\n"))
# splitting into two groups
TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]] <- cut2(TEMP.DF[,TRAITS.TARGET.RANK[target_of_interest]], g = 2)
cat(paste0(" > cross tabulation of ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
show(table(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]))
cat(paste0("\n > fitting the model for ",TRAITS.TARGET.RANK[target_of_interest],"-stratum.\n"))
fit <- survfit(as.formula(paste0("y ~ ", TRAITS.TARGET.RANK[target_of_interest])), data = TEMP.DF)
cat(paste0("\n > make a Kaplan-Meier-shizzle...\n"))
# make Kaplan-Meier curve and save it
show(ggsurvplot(fit, data = TEMP.DF,
palette = c("#DB003F", "#1290D9"),
# palete = c("F59D10", "#DB003F", "#49A01D", "#1290D9"),
linetype = c(1,2),
ylim = c(0.75, 1),
# linetype = c(1,2,3,4),
# conf.int = FALSE, conf.int.fill = "#595A5C", conf.int.alpha = 0.1,
pval = FALSE, pval.method = FALSE, pval.size = 4,
risk.table = TRUE, risk.table.y.text = FALSE, tables.y.text.col = TRUE, fontsize = 4,
censor = FALSE,
legend = "right",
legend.title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
legend.labs = c("low", "high"),
title = paste0("Risk of ",ep,""), xlab = "Time [days]", font.main = c(16, "bold", "black")))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.survival.",ep,".2G.",
TRAITS.TARGET.RANK[target_of_interest],".90days.pdf"), width = 12, height = 10, onefile = FALSE)
cat(paste0("\n > perform the Cox-regression fashizzle and plot it...\n"))
### Do Cox-regression and plot it
### MODEL 2 adjusted for age, sex, hypertension, diabetes, smoking, LDL-C levels, lipid-lowering drugs, antiplatelet drugs, eGFR, BMI, history of CVD, level of stenosis
cox = coxph(Surv(TEMP.DF[,eptime], event) ~ TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]]+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
coxplot = coxph(Surv(TEMP.DF[,eptime], event) ~ strata(TEMP.DF[[ TRAITS.TARGET.RANK[target_of_interest] ]])+Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose, data = TEMP.DF)
plot(survfit(coxplot), main = paste0("Cox proportional hazard of [",ep,"] per [",eptime,"]."),
ylim = c(0.75, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
# ylim = c(0, 1), xlim = c(0,3), col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
ylab = "Suvival probability", xlab = "FU time [days]",
mark.time = FALSE, axes = FALSE, bty = "n")
legend("topright",
c("low", "high"),
title = paste0("",TRAITS.TARGET.RANK[target_of_interest],""),
col = c("#DB003F", "#1290D9"),
lty = c(1,2), lwd = 2,
bty = "n")
axis(side = 1, at = seq(0, 3, by = 1))
axis(side = 2, at = seq(0, 1, by = 0.2))
dev.copy2pdf(file = paste0(COX_loc,"/",
Today,".AERNASE.clin.hdac9.Cox.",ep,".2G.",
# Today,".AERNASE.clin.hdac9.Cox.",ep,".4G.",
TRAITS.TARGET.RANK[target_of_interest],".MODEL2.90days.pdf"), height = 12, width = 10, onefile = TRUE)
show(summary(cox))
cat(paste0("\n > writing the Cox-regression fashizzle to Excel...\n"))
COX.results.TEMP <- data.frame(matrix(NA, ncol = 12, nrow = 0))
COX.results.TEMP[1,] = COX.STAT(cox, "AERNASE.clin.hdac9", ep, TRAITS.TARGET.RANK[target_of_interest])
COX.results = rbind(COX.results, COX.results.TEMP)
}
}
cat("- Edit the column names...\n")
colnames(COX.results) = c("Dataset", "Outcome", "CpG",
"Beta", "s.e.m.",
"HR", "low95CI", "up95CI",
"Z-value", "P-value", "SampleSize", "N_events")
cat("- Correct the variable types...\n")
COX.results$Beta <- as.numeric(COX.results$Beta)
COX.results$s.e.m. <- as.numeric(COX.results$s.e.m.)
COX.results$HR <- as.numeric(COX.results$HR)
COX.results$low95CI <- as.numeric(COX.results$low95CI)
COX.results$up95CI <- as.numeric(COX.results$up95CI)
COX.results$`Z-value` <- as.numeric(COX.results$`Z-value`)
COX.results$`P-value` <- as.numeric(COX.results$`P-value`)
COX.results$SampleSize <- as.numeric(COX.results$SampleSize)
COX.results$N_events <- as.numeric(COX.results$N_events)
AERNASE.clin.hdac9.COX.results <- COX.results
# Save the data
cat("- Writing results to Excel-file...\n")
head.style <- createStyle(textDecoration = "BOLD")
write.xlsx(AERNASE.clin.hdac9.COX.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Cox.2G.MODEL2.90days.xlsx"),
creator = "Sander W. van der Laan",
sheetName = "Results", headerStyle = head.style,
rowNames = FALSE, colNames = TRUE, overwrite = TRUE)
# Removing intermediates
cat("- Removing intermediate files...\n")
rm(TEMP.DF, target_of_interest, fit, cox, coxplot, COX.results, COX.results.TEMP, head.style, AERNASE.clin.hdac9.COX.results)We correlated plaque levels of the biomarkers.
# Installation of ggcorrplot()
# --------------------------------
# if(!require(devtools))
# install.packages("devtools")
# devtools::install_github("kassambara/ggcorrplot")
library(ggcorrplot)
# Creating matrix - inverse-rank transformation
# --------------------------------
# AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9,
# select = c("IL6_rank", "MCP1_rank", "IL6_pg_ug_2015_rank", "MCP1_pg_ug_2015_rank", "IL6R_pg_ug_2015_rank",
# TRAITS.BIN, TRAITS.CON.RANK)
# )
# AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9,
# select = c("MCP1_rank", "MCP1_pg_ug_2015_rank",
# TRAITS.BIN, TRAITS.CON.RANK)
# )
AERNASE.clin.hdac9.temp <- subset(AERNASE.clin.hdac9,
select = c("HDAC9",
TRAITS.BIN,
TRAITS.CON.RANK,
"Symptoms.5G", "AsymptSympt", "EP_major", "EP_composite")
)
AERNASE.clin.hdac9.temp$CalcificationPlaque <- as.numeric(AERNASE.clin.hdac9.temp$CalcificationPlaque)
AERNASE.clin.hdac9.temp$CollagenPlaque <- as.numeric(AERNASE.clin.hdac9.temp$CollagenPlaque)
AERNASE.clin.hdac9.temp$Fat10Perc <- as.numeric(AERNASE.clin.hdac9.temp$Fat10Perc)
AERNASE.clin.hdac9.temp$MAC_binned <- as.numeric(AERNASE.clin.hdac9.temp$MAC_binned)
AERNASE.clin.hdac9.temp$SMC_binned <- as.numeric(AERNASE.clin.hdac9.temp$SMC_binned)
AERNASE.clin.hdac9.temp$IPH <- as.numeric(AERNASE.clin.hdac9.temp$IPH)
AERNASE.clin.hdac9.temp$Symptoms.5G <- as.numeric(AERNASE.clin.hdac9.temp$Symptoms.5G)
AERNASE.clin.hdac9.temp$AsymptSympt <- as.numeric(AERNASE.clin.hdac9.temp$AsymptSympt)
AERNASE.clin.hdac9.temp$EP_major <- as.numeric(AERNASE.clin.hdac9.temp$EP_major)
AERNASE.clin.hdac9.temp$EP_composite <- as.numeric(AERNASE.clin.hdac9.temp$EP_composite)
str(AERNASE.clin.hdac9.temp)'data.frame': 623 obs. of 15 variables:
$ HDAC9 : int 24 11 7 5 13 46 14 0 30 12 ...
$ CalcificationPlaque : num 2 2 2 1 1 1 2 2 2 1 ...
$ CollagenPlaque : num 2 2 2 2 2 2 1 2 2 2 ...
$ Fat10Perc : num 2 2 2 2 1 1 2 1 2 1 ...
$ IPH : num 1 2 2 1 1 1 1 2 1 1 ...
$ MAC_binned : num 2 1 2 1 1 2 2 1 2 1 ...
$ SMC_binned : num 2 2 2 2 2 2 1 2 2 2 ...
$ MAC_rankNorm : num -0.0695 -0.7504 0.9099 -0.3812 -1.0959 ...
$ SMC_rankNorm : num 1.526 0.58 0.674 1.26 0.341 ...
$ MAC_SMC_ratio_rank : num -0.215 -0.888 0.722 -0.366 -1.331 ...
$ VesselDensity_rankNorm: num -0.723 0.769 0.854 1.14 -0.861 ...
$ Symptoms.5G : num 6 6 5 6 2 5 1 1 1 1 ...
$ AsymptSympt : num 3 3 3 3 2 3 1 1 1 1 ...
$ EP_major : num 2 2 2 2 3 3 3 2 2 2 ...
$ EP_composite : num 2 3 3 2 3 3 3 2 2 2 ...AERNASE.clin.hdac9.matrix.RANK <- as.matrix(AERNASE.clin.hdac9.temp)
rm(AERNASE.clin.hdac9.temp)corr_biomarkers.rank <- round(cor(AERNASE.clin.hdac9.matrix.RANK,
use = "pairwise.complete.obs", #the correlation or covariance between each pair of variables is computed using all complete pairs of observations on those variables
method = "spearman"), 3)
corr_biomarkers.rank HDAC9 CalcificationPlaque CollagenPlaque Fat10Perc IPH MAC_binned SMC_binned MAC_rankNorm SMC_rankNorm MAC_SMC_ratio_rank
HDAC9 1.000 0.045 -0.069 0.109 0.055 0.045 -0.091 0.000 -0.140 -0.084
CalcificationPlaque 0.045 1.000 0.126 0.029 0.083 -0.079 0.068 -0.006 0.015 0.041
CollagenPlaque -0.069 0.126 1.000 -0.077 0.006 0.008 0.334 0.015 0.303 0.091
Fat10Perc 0.109 0.029 -0.077 1.000 0.201 0.293 -0.154 0.197 -0.250 -0.012
IPH 0.055 0.083 0.006 0.201 1.000 0.088 -0.120 0.128 -0.080 -0.033
MAC_binned 0.045 -0.079 0.008 0.293 0.088 1.000 0.042 0.627 -0.041 -0.046
SMC_binned -0.091 0.068 0.334 -0.154 -0.120 0.042 1.000 0.112 0.557 0.089
MAC_rankNorm 0.000 -0.006 0.015 0.197 0.128 0.627 0.112 1.000 0.212 -0.015
SMC_rankNorm -0.140 0.015 0.303 -0.250 -0.080 -0.041 0.557 0.212 1.000 0.182
MAC_SMC_ratio_rank -0.084 0.041 0.091 -0.012 -0.033 -0.046 0.089 -0.015 0.182 1.000
VesselDensity_rankNorm -0.034 -0.050 0.171 0.014 0.203 0.073 0.109 0.085 0.077 0.026
Symptoms.5G 0.069 -0.098 -0.107 0.139 0.022 0.065 -0.126 0.039 -0.172 -0.002
AsymptSympt 0.074 -0.113 -0.101 0.144 0.031 0.067 -0.172 0.023 -0.186 -0.007
EP_major -0.058 0.062 0.047 -0.017 0.004 -0.041 0.078 -0.003 0.034 0.010
EP_composite -0.008 0.066 0.049 -0.036 0.075 -0.063 0.031 -0.005 0.042 0.047
VesselDensity_rankNorm Symptoms.5G AsymptSympt EP_major EP_composite
HDAC9 -0.034 0.069 0.074 -0.058 -0.008
CalcificationPlaque -0.050 -0.098 -0.113 0.062 0.066
CollagenPlaque 0.171 -0.107 -0.101 0.047 0.049
Fat10Perc 0.014 0.139 0.144 -0.017 -0.036
IPH 0.203 0.022 0.031 0.004 0.075
MAC_binned 0.073 0.065 0.067 -0.041 -0.063
SMC_binned 0.109 -0.126 -0.172 0.078 0.031
MAC_rankNorm 0.085 0.039 0.023 -0.003 -0.005
SMC_rankNorm 0.077 -0.172 -0.186 0.034 0.042
MAC_SMC_ratio_rank 0.026 -0.002 -0.007 0.010 0.047
VesselDensity_rankNorm 1.000 -0.042 -0.042 0.111 0.167
Symptoms.5G -0.042 1.000 0.885 0.004 -0.047
AsymptSympt -0.042 0.885 1.000 0.002 -0.051
EP_major 0.111 0.004 0.002 1.000 0.675
EP_composite 0.167 -0.047 -0.051 0.675 1.000corr_biomarkers_p.rank <- ggcorrplot::cor_pmat(AERNASE.clin.hdac9.matrix.RANK, use = "pairwise.complete.obs", method = "spearman", exact = FALSE)# Add correlation coefficients
# --------------------------------
# argument lab = TRUE
ggcorrplot(corr_biomarkers.rank,
method = "square",
type = "lower",
title = "Cross biomarker correlations",
show.legend = TRUE, legend.title = bquote("Spearman's"~italic(rho)),
ggtheme = ggplot2::theme_minimal, outline.color = "#FFFFFF",
show.diag = TRUE,
hc.order = FALSE,
lab = FALSE,
digits = 3,
# p.mat = corr_biomarkers_p.rank, sig.level = 0.05,
colors = c("#1290D9", "#FFFFFF", "#E55738"))
library(data.table)
# flattenCorrMatrix
# --------------------------------
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
biomarker_row = rownames(cormat)[row(cormat)[ut]],
biomarker_column = rownames(cormat)[col(cormat)[ut]],
spearman_cor =(cormat)[ut],
pval = pmat[ut]
)
}
corr_biomarkers.rank.df <- as.data.table(flattenCorrMatrix(corr_biomarkers.rank, corr_biomarkers_p.rank))
DT::datatable(corr_biomarkers.rank.df)NA# chart of a correlation matrix
# --------------------------------
# Alternative solution https://www.r-graph-gallery.com/199-correlation-matrix-with-ggally.html
install.packages.auto("PerformanceAnalytics")
chart.Correlation.new <- function (R, histogram = TRUE, method = c("pearson", "kendall",
"spearman"), ...)
{
x = checkData(R, method = "matrix")
if (missing(method))
method = method[1]
cormeth <- method
panel.cor <- function(x, y, digits = 2, prefix = "", use = "pairwise.complete.obs",
method = cormeth, cex.cor, ...) {
usr <- par("usr")
on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- cor(x, y, use = use, method = method)
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste(prefix, txt, sep = "")
if (missing(cex.cor))
cex <- 0.8/strwidth(txt)
test <- cor.test(as.numeric(x), as.numeric(y), method = method)
Signif <- symnum(test$p.value, corr = FALSE, na = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), symbols = c("***",
"**", "*", ".", " "))
text(0.5, 0.5, txt, cex = cex * (abs(r) + 0.3)/1.3)
text(0.8, 0.8, Signif, cex = cex, col = 2)
}
f <- function(t) {
dnorm(t, mean = mean(x), sd = sd.xts(x))
}
dotargs <- list(...)
dotargs$method <- NULL
rm(method)
hist.panel = function(x, ... = NULL) {
par(new = TRUE)
hist(x, col = "#1290D9", probability = TRUE, axes = FALSE,
# hist(x, col = "light gray", probability = TRUE, axes = FALSE,
main = "", breaks = "FD")
lines(density(x, na.rm = TRUE), col = "#E55738", lwd = 1)
rug(x)
}
if (histogram)
pairs(x, gap = 0, lower.panel = panel.smooth, upper.panel = panel.cor,
diag.panel = hist.panel, ...)
else pairs(x, gap = 0, lower.panel = panel.smooth, upper.panel = panel.cor, ...)
}
chart.Correlation.new(AERNASE.clin.hdac9.matrix.RANK, method = "spearman", histogram = TRUE, pch = 3)
# alternative chart of a correlation matrix
# --------------------------------
# Alternative solution https://www.r-graph-gallery.com/199-correlation-matrix-with-ggally.html
install.packages.auto("GGally")
# Quick display of two cabapilities of GGally, to assess the distribution and correlation of variables
library(GGally)
# From the help page:
ggpairs(AERNASE.clin.hdac9,
columns = c("HDAC9", TRAITS.BIN, TRAITS.CON.RANK, "Symptoms.5G", "AsymptSympt", "EP_major", "EP_composite"),
columnLabels = c("HDAC9",
"Calcification", "Collagen", "Fat 10%", "IPH", "Macrophages (binned)", "SMC (binned)", "Macrophages", "SMC", "Macrophage/SMC", "Vessel density",
"Symptoms", "Symptoms (grouped)", "MACE", "Composite"),
method = c("spearman"),
# ggplot2::aes(colour = Gender),
progress = FALSE)`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
For the purpose of downstream analyses we save the
AEDB.CEA object, and the HDAC9 objects.
saveRDS(AEDB.CEA, file = paste0(OUT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".AEDB.CEA.withupdates.RDS"))
saveRDS(AEDB.CEA, file = paste0(OUT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".AERNASE.clin.hdac9.df.RDS"))
saveRDS(AEDB.CEA, file = paste0(OUT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".AERNASE.hdac9.se.RDS"))
saveRDS(AEDB.CEA, file = paste0(OUT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".AERNASE.hdac9.genedata.RDS"))Version: v1.0.1
Last update: 2021-03-19
Written by: Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description: Script to analyse HDAC9 from the Ather-Express Biobank Study.
Minimum requirements: R version 3.5.2 (2018-12-20) -- 'Eggshell Igloo', macOS Mojave (10.14.2).
**MoSCoW To-Do List**
The things we Must, Should, Could, and Would have given the time we have.
_M_
_S_
_C_
_W_
**Changes log**
* v1.0.1 Update to main AEDB (there is an error in the Age-variable in the new version). More patients in the bulk (623 vs 611 with the newer dataset).
* v1.0.0 Inital version.sessionInfo()R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Monterey 12.2.1
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats4 grid tools stats graphics grDevices utils datasets methods base
other attached packages:
[1] RColorBrewer_1.1-2 SeuratObject_4.0.4 Seurat_4.0.6
[4] PerformanceAnalytics_2.0.4 xts_0.12.1 zoo_1.8-9
[7] Hmisc_4.6-0 Formula_1.2-4 lattice_0.20-45
[10] survminer_0.4.9 survival_3.3-1 MASS_7.3-54
[13] ggcorrplot_0.1.3.999 GGally_2.1.2 annotables_0.1.91
[16] EnhancedVolcano_1.12.0 ggrepel_0.9.1 AnnotationFilter_1.18.0
[19] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 mygene_1.30.0 org.Hs.eg.db_3.14.0
[22] DESeq2_1.34.0 SummarizedExperiment_1.24.0 MatrixGenerics_1.6.0
[25] matrixStats_0.61.0 GenomicFeatures_1.46.3 AnnotationDbi_1.56.2
[28] Biobase_2.54.0 GenomicRanges_1.46.1 GenomeInfoDb_1.30.0
[31] IRanges_2.28.0 S4Vectors_0.32.3 BiocGenerics_0.40.0
[34] patchwork_1.1.0.9000 labelled_2.9.0 openxlsx_4.2.5
[37] sjPlot_2.8.10 UpSetR_1.4.0 ggpubr_0.4.0
[40] forestplot_2.0.1 checkmate_2.0.0 magrittr_2.0.2
[43] pheatmap_1.0.12 devtools_2.4.3 usethis_2.1.5
[46] BlandAltmanLeh_0.3.1 tableone_0.13.0 haven_2.4.3
[49] eeptools_1.2.4 DT_0.20 knitr_1.37
[52] forcats_0.5.1 stringr_1.4.0 purrr_0.3.4
[55] tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1
[58] data.table_1.14.2 naniar_0.6.1 tidyr_1.1.4
[61] dplyr_1.0.7 optparse_1.7.1 readr_2.1.1
[64] pander_0.6.4 rmarkdown_2.11 worcs_0.1.9.1
loaded via a namespace (and not attached):
[1] mitools_2.4 pbapply_1.5-0 vctrs_0.3.8 mgcv_1.8-39 blob_1.2.2 spatstat.data_2.1-2
[7] later_1.3.0 nloptr_1.2.2.3 DBI_1.1.2 uwot_0.1.11 rappdirs_0.3.3 gsubfn_0.7
[13] jpeg_0.1-9 zlibbioc_1.40.0 sjmisc_2.8.9 htmlwidgets_1.5.4 mvtnorm_1.1-3 future_1.23.0
[19] leiden_0.3.9 parallel_4.1.2 irlba_2.3.5 markdown_1.1 Rcpp_1.0.8.3 KernSmooth_2.23-20
[25] promises_1.2.0.1 limma_3.50.0 DelayedArray_0.20.0 ggeffects_1.1.1 pkgload_1.2.4 fs_1.5.2
[31] textshaping_0.3.6 ranger_0.13.1 digest_0.6.29 png_0.1-7 sctransform_0.3.2 cowplot_1.1.1
[37] pkgconfig_2.0.3 ggbeeswarm_0.6.0 estimability_1.3 minqa_1.2.4 reticulate_1.22 beeswarm_0.4.0
[43] xfun_0.29 bslib_0.3.1 tidyselect_1.1.1 performance_0.8.0 reshape2_1.4.4 ica_1.0-2
[49] viridisLite_0.4.0 rtracklayer_1.54.0 pkgbuild_1.3.1 rlang_1.0.2 jquerylib_0.1.4 glue_1.6.2
[55] modelr_0.1.8 emmeans_1.7.2 ggsignif_0.6.3 bayestestR_0.11.5 labeling_0.4.2 maptools_1.1-2
[61] httpuv_1.6.5 class_7.3-20 Rttf2pt1_1.3.9 TH.data_1.1-0 annotate_1.72.0 jsonlite_1.7.2
[67] XVector_0.34.0 bit_4.0.4 mime_0.12 systemfonts_1.0.3 gridExtra_2.3 Rsamtools_2.10.0
[73] stringi_1.7.6 insight_0.16.0 processx_3.5.2 spatstat.sparse_2.1-0 scattermore_0.7 survey_4.1-1
[79] quadprog_1.5-8 bitops_1.0-7 cli_3.2.0 sqldf_0.4-11 maps_3.4.0 RSQLite_2.2.9
[85] prereg_0.5.0 rticles_0.22 rsconnect_0.8.25 rstudioapi_0.13 GenomicAlignments_1.30.0 nlme_3.1-155
[91] locfit_1.5-9.4 listenv_0.8.0 miniUI_0.1.1.1 survMisc_0.5.5 dbplyr_2.1.1 sessioninfo_1.2.2
[97] readxl_1.3.1 lifecycle_1.0.1 munsell_0.5.0 cellranger_1.1.0 ggsci_2.9 codetools_0.2-18
[103] coda_0.19-4 vipor_0.4.5 lmtest_0.9-39 sys_3.4 htmlTable_2.4.0 proto_1.0.0
[109] xtable_1.8-4 ROCR_1.0-11 abind_1.4-5 farver_2.1.0 parallelly_1.30.0 km.ci_0.5-2
[115] credentials_1.3.2 RANN_2.6.1 askpass_1.1 visdat_0.5.3 BiocIO_1.4.0 sjstats_0.18.1
[121] goftest_1.2-3 RcppAnnoy_0.0.19 cluster_2.1.2 future.apply_1.8.1 extrafontdb_1.0 Matrix_1.4-0
[127] ellipsis_0.3.2 prettyunits_1.1.1 lubridate_1.8.0 ggridges_0.5.3 reprex_2.0.1 igraph_1.2.11
[133] sjlabelled_1.1.8 remotes_2.4.2 parameters_0.17.0 spatstat.utils_2.3-0 testthat_3.1.1 getopt_1.20.3
[139] htmltools_0.5.2 BiocFileCache_2.2.0 yaml_2.2.1 utf8_1.2.2 plotly_4.10.0 XML_3.99-0.8
[145] e1071_1.7-9 foreign_0.8-82 withr_2.5.0 fitdistrplus_1.1-6 BiocParallel_1.28.3 bit64_4.0.5
[151] effectsize_0.6.0.1 multcomp_1.4-18 ProtGenerics_1.26.0 spatstat.core_2.3-2 Biostrings_2.62.0 ragg_1.2.1
[157] memoise_2.0.1 evaluate_0.14 geneplotter_1.72.0 tzdb_0.2.0 extrafont_0.17 callr_3.7.0
[163] ps_1.6.0 curl_4.3.2 fansi_1.0.2 tensor_1.5 cachem_1.0.6 desc_1.4.0
[169] deldir_1.0-6 proj4_1.0-10.1 rjson_0.2.21 rstatix_0.7.0 rprojroot_2.0.2 sass_0.4.0
[175] sandwich_3.0-1 RCurl_1.98-1.5 proxy_0.4-26 car_3.0-12 xml2_1.3.3 httr_1.4.2
[181] assertthat_0.2.1 boot_1.3-28 globals_0.14.0 R6_2.5.1 nnet_7.3-17 progress_1.2.2
[187] genefilter_1.76.0 KEGGREST_1.34.0 ggrastr_1.0.1 splines_4.1.2 carData_3.0-5 colorspace_2.0-3
[193] generics_0.1.1 base64enc_0.1-3 chron_2.3-56 gridtext_0.1.4 gert_1.5.0 pillar_1.7.0
[199] ggalt_0.4.0 sp_1.4-6 GenomeInfoDbData_1.2.7 plyr_1.8.6 gtable_0.3.0 rvest_1.0.2
[205] zip_2.2.0 restfulr_0.0.13 latticeExtra_0.6-29 biomaRt_2.50.2 fastmap_1.1.0 Cairo_1.5-14
[211] crosstalk_1.2.0 datawizard_0.3.0 vcd_1.4-9 broom_0.7.11 openssl_1.4.6 scales_1.1.1
[217] arm_1.12-2 filelock_1.0.2 backports_1.4.1 lme4_1.1-27.1 hms_1.1.1 Rtsne_0.15
[223] shiny_1.7.1 KMsurv_0.1-5 ash_1.0-15 polyclip_1.10-0 lazyeval_0.2.2 crayon_1.5.0
[229] reshape_0.8.8 rpart_4.1.16 spatstat.geom_2.3-1 compiler_4.1.2 ggtext_0.1.1 save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".bulkRNAseq.main_analysis.RData"))| © 1979-2022 Sander W. van der Laan | s.w.vanderlaan[at]gmail.com | swvanderlaan.github.io. |